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Abstract 

We introduce a new efficient algorithm for phase space generation. A parton shower is used to 
distribute events across all of multiplicity, flavor, and phase space, and these events can then 
be reweighted to any desired analytic distribution. To verify this method, we reproduce the 
Q-f e + e~ — ► n jets tree-level result of traditional matrix element tools. We also show how to improve 

tree-level matrix elements automatically with leading-logarithmic resummation. This algorithm 
is particularly useful in the context of a new framework for event generation called GenEvA. In 
a companion paper [Q], we show how the GenEvA framework can address contemporary issues in 
' <— ] , event generation. 
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1. INTRODUCTION 



A. Motivation 



Precision theoretical calculations are vital for any comparison between experimental mea- 
surements and underlying theoretical models. While in many cases we are interested in the 
effects of physics Beyond the Standard Model, a precise understanding of Standard Model 
backgrounds is always required for such a comparison to be useful. To be able to study 
the effect of experimental cuts and detector effects, the most useful form for theoretical 
predictions is as fully exclusive events with a high multiplicity of hadrons in the final state. 

Ignoring the issue of hadronization and other nonperturbative effects, the challenge for 
constructing exclusive event generators at the perturbative level is that full quantum field 
theoretic results for high-multiplicity final states become prohibitively difficult both to cal- 
culate and to implement. At tree level, the number of required Feynman diagrams grows 
factorially with the number of partons in the final state, so that a direct calculation is not 
feasible. Even if an expression for the partonic calculation can be obtained, the large volume 
of high-multiplicity phase space makes it difficult to efficiently distribute events according 
to those distributions. Both of these constraints imply that in practice, it is impossible 
to use full quantum information to calculate typical partonic final states arising in collider 
environments. 

It has been known for some time that for emissions at small transverse momentum, much 
of the quantum nature of QCD becomes subdominant, such that emissions in QCD can 
be approximated by classical splitting functions @, [| ^j. These splitting functions form 
the basis for so-called parton showers, which build up arbitrarily complicated final states 
by starting with low-multiplicity final states and recursively splitting one final state parton 
into two. The recursive nature of parton showers allows them to be cast into a Markov 
chain algorithm. Parton showers form the basis for multi-purpose event generators like 
Pythia f| ||, Herwig 0, g, g], Ariadne 0, Isajet fy], and Sherpa fT^j, which build up 
QCD-like events from simple hard-scattering matrix elements. 

While parton shower algorithms easily generate high-multiplicity events, the resulting 
distributions are not based on full QCD, and for example do not accurately describe emis- 
sions at large transverse momentum. For this reason event generators based on full ma- 
trix element information have been constructed. Examples for generators based on tree- 
level calculations are Alpgen @, MadEvent Q, CompHEP |L5|, AMEGIC |L6f, Whizard |L7 



and Helac-Phegas |J, fL9|, 131, wnile 

programs using next-to-leading order (NLO) calcu- 
lations also exist 



21 , NLO Jet 22 , PHOX 23 , and VBFNLO 



In 



examples being MCFM 

the last few years, significant progress has been made towards combining full QCD cal 
culations, distributed with traditional phase space generators, with parton shower algo- 

EEL E71 ra, m 



rithms 25 



31], B2, |33|, p4l |35|, B6l |37| . Encouraged by these successes, it is 



worth considering whether additional progress could be made by considering a fundamentally 
new approach to event generation. 

In this paper, we introduce a novel use of the parton shower as a phase space generator 
for generic partonic calculations. The idea is to generate an event with a given final state 
multiplicity using a simple parton shower algorithm, and then to reweight this event to any 
desired partonic distribution. This algorithm is a key component of the GenEvA framework 
for event generation that we discuss in a companion paper Pj] , so we will refer to the phase 
space generator presented here as the GenEvA algorithm — for Generate Events Analytically. 
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The GenEvA framework asserts that the choice of partonic calculations should not be 
determined by algorithmic needs but rather by the underlying full theory information we 
want to describe. For this to be possible, the GenEvA framework requires a "projectable" 
phase space generator with a built-in notion of a matching scale, and the GenEvA algorithm 
is an example of such a generator. In the companion paper [|IJ], we implement the following 
partonic results for e + e~ — > n jets: 

1. LO/LL: Tree-level (LO) matrix elements with up to 6 final state partons, improved 
by leading logarithms (LL), 

2. NLO/LL: Next-to- leading order (NLO) matrix elements including the virtual 2- 
parton diagram, improved by leading logarithms, 

3. NLO/LO/LL: A combination of the above two results, where the NLO calcula- 
tions are supplemented by higher-multiplicity tree-level calculation, all with leading- 
logarithmic resummation. 



While both LO/LL [EE EEL E7L Eg EM I3TL B3, B3, B§ m M, O and NLO/LL [BR m ETC 



§T], 0, H| H EH E§ E3> EH © H EO results nave been implemented elsewhere, the third 
result has never been implemented before. Furthermore, these distributions are obtained 
without the need for negative-weight events. 

From a technical point of view, reweighting a parton shower has become possible with the 



introduction of the analytic parton shower in Ref. ||52|| . There, the probabilities to obtain 
individual branches are independent, so that the total probability of generating a given 
shower history is just the product of multiple splitting probabilities. This allows one to 
analytically calculate the exact probability for generating a given event. Because the parton 
shower covers Lorentz-invariant phase space multiple times, we also introduce a method to 
eliminate this double-counting problem while making maximal use of the full event sample 
given by the parton shower. 

Some of the techniques in this paper are readily applicable to existing Monte Carlo pro- 
grams, and could be implemented via a simple change in bookkeeping instead of a complete 
code rewrite. Indeed, the idea of reweighting parton showers to matrix elements appeared 
in the literature over 20 years ago |26], |53|, [53], F55], [56], |57|, |58|, [)9|] and is implemented in 



Monte Carlo programs like Pythia |j], |6[]. The GenEvA algorithm offers a generalization of 
this well-known method. As we will see, some algorithmic and conceptual ideas are taken 
directly from ALPGEN |I3| and MadEvent |T4[], so the code optimizations present in more 



mature Monte Carlo programs could aid in further increasing the efficiency of the GenEvA 
algorithm. 

As in the companion paper jTJ, we will focus on the process e + e~ — *> n jets, though in the 
conclusions in Sec. |] we discuss obvious generalizations to hadronic collisions with heavy 
resonances. In Sec. we introduce the main ingredients of the GenEvA algorithm: the parton 
shower probability, the Jacobian transformation to Lorentz-invariant phase space, the over- 
counting factor, and the desired matrix element. Sees. f|, |5], and || discuss these ingredients 
in turn, and Sec. |] also contains a discussion how the structure of the GenEvA algorithm 
allows for an automatic leading-logarithmic improvement of tree-level matrix elements. In 
Sec. |7| we verify that the GenEvA algorithm works as a proper phase space generator, showing 
as an example that it can reproduce the results of traditional tree-level event generators, 
and discuss the efficiency of the algorithm. The reader not interested in the technical details 
may skip Sees. |^ through || and directly go from Sec. to Sec. |7|. 
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Before going into the details of the GenEvA algorithm, we first want to discuss why the 
parton shower is an appropriate and beneficial starting point for event generation. 



B. Phase Space from a Parton Shower 

What does it mean to generate phase space from a parton shower? Ignoring the fact 
that the parton shower is derived from the soft-collinear limit of QCD, there are three 
main distinguishing features. First, it is well known that n-body phase space can be built 
recursively from (n — l)-body phase space, and a parton shower is an explicit implementation 
of this recursion. Second, a single parton shower run covers all of multiplicity, flavor, and 
phase space (hereafter referred to just as "phase space") in one shot, which is an example 
of "single-channel" integration. Third, the parton shower preserves probability, meaning 
that during an event run, the fraction (or "grid") of events of a given multiplicity, flavor, 
and kinematic structure will remain constant. Thus, as a phase space generator, the parton 
shower is a recursive, single-channel, fixed-grid algorithm. 

The GenEvA algorithm is then diametrically opposed to modern phase space integra- 
tors like MadEvent [H|, which use non-recursive, multi-channel, adaptive-grid algorithms, 



based on Vegas pO] and the Bases/Spring package pM. While event generators like 



CalcHEP/CompHEP ^2] do use recursive methods (though not a parton shower) to gen- 
erate phase space, algorithms such as RAMBO |)3| and SARGE |)4j were developed in order to 
avoid using a recursive definition of phase space because of problems of efficiency Almost 
every phase space generator uses some notion of adaptation, in order to avoid the problem of 
large over- or under-sampling of phase space regions. And advanced phase space algorithms 
use multi-channel integration to create a single event sample culled from separate runs. 

To see why GenEvA is a step forward rather than a step backward, it is important to 
note one crucial commonality between MadEvent and GenEvA: Both know about the singu- 
larity structure of QCD from 1/Q 2 poles. In practice, we will see that for e + e~ — > n jets, 
GenEvA performs at or above MadEvent efficiencies. This suggests that the use of advanced 
non- recursive, multi-channel, adaptive-grid algorithms by other phase space generators is 
compensating for something that GenEvA has already built-in, namely, an excellent approx- 
imation to QCD-like events through the parton shower. In other words, despite the fact 
that the recursive, single-channel, fixed-grid machinery of GenEvA is so primitive, the par- 
ton shower is sufficiently close to the right answer that more advanced methods are simply 
not necessary. Whether this fact will continue to be true in hadronic collisions is an open 
question. 

While the machinery of GenEvA is simple-minded, the technical details needed to turn 
a parton shower into a phase space generator are not. The fact that the parton shower 
preserves probability means that one needs a thorough understanding of exactly what that 
probability distribution is, and the bulk of this paper is devoted to explaining in detail how 



the analytic parton shower in Ref. |52[ can become an efficient generator of Lorentz-invariant 
phase space. 

For example, a crucial technical detail is that since the parton shower generates the entire 
phase space out of successive 1 — >• 2 splittings, it will usually distribute events with more 
partons in the final state than we have full QCD calculations available for. The GenEvA 
algorithm allows such events to be "truncated" to lower-multiplicity final states in a way 
that still allows for an analytic calculation of the resulting truncation probability. This is a 
concrete realization of the "phase space projection" requirement of the GenEvA framework M, 
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though the form implemented here is a "reversible projection" that has additional analytic 
properties beyond what is strictly speaking necessary in that context. 



C. Benefits of GenEvA 

There are numerous benefits to using the parton shower as a phase space generator, and 
many of these benefits actually derive from the fact that GenEvA is a recursive, single-channel, 
fixed-grid algorithm. 

First, the parton shower is based on the symmetry structures of QCD. All possible flavor 
and color structures are generated with approximately the right probabilities, meaning that 
processes that are naturally linked, such as e + e~ — > qq + n gluons and e + e~ — > qqq'q' + 
m gluons, need not have separate phase space integrations. The single-channel nature of 
the GenEvA algorithm makes the program trivial to parallelize, because a single GenEvA run 
contains a complete sampling of all desired physics processes. 

Second, the parton shower is based on the singularity structures of QCD. GenEvA repro- 
duces both the singular pole and leading-logarithmic behavior of full QCD distributions, 
yielding events with moderate weights. More importantly, the program distributes events 
with appropriate Sudakov factors already, allowing for an automatic leading-log im- 



provement of any tree-level matrix element. Of course, these generated Sudakov factors can 
also be calculated analytically and removed or exchanged if desired. 

Third, because the phase space generation is recursive, a single event can have a nested 
set of consistent n-body phase space variables. This allows considerable flexibility in the 
types of parton-level calculations that can be handled, and allows a simple implementation 
of the NLO/LO/LL merged sample from the companion paper 

Fourth, a single event can have numerous different weights corresponding to different 
theoretical assumptions and tunable parameters. Using the feature of truncation, GenEvA 
allows the same event to be interpreted as coming from an n-body or m-body matrix element 
with n not necessarily the same as m. 1 This makes it possible to have a single set of events 
that is capable of describing multiple theoretical distributions, greatly simplifying the process 
of parton-level systematic error estimation. 

Moreover, by running this single event set through the required detector simulation, one 
obtains the equivalent of multiple detector-improved distributions, each with roughly the 
statistical significance of the full number of events. Given the large computational cost of 
full detector simulation, it is highly desirable to reuse Monte Carlo in this way. 2 Because 
the GenEvA algorithm includes the proper singularities and Sudakov factors out-of-the-box, 
this reweighting strategy will potentially be feasible even on large Monte Carlo sets. 

Finally, GenEvA is well-suited to situations where it is difficult or tedious to distribute 
events according to a certain kind of theoretical calculation. For example, the application 
of soft-collinear effective theory (SCET) |j6], |67|, |6{| to event generation j7(| fn], |72| is 



most useful once the final state kinematics and matching scales fx are already known. If we 
are interested in calculating subleading logarithms using SCET, it is hard to analytically 



1 In a traditional phase space generator, it is possible to reweight events distributed according to one n-body 
matrix element to a different n-body matrix element, but the value of n must be fixed. 

2 Because detector simulation is a linear operation, it is even possible to try different theoretical distributions 
after detector simulation. 
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determine the proper SCET amplitudes over all of phase space. However, as a template 
for QCD-like events, GenEvA determines the kinematics and scales without reference to the 
amplitudes, allowing in principle numeric methods to be used to calculate SCET operator 
mixing. 

Because of the above benefits, it is tempting to identify the GenEvA framework introduced 
in the companion paper [|TJ with the GenEvA algorithm that will be presented here. This 
is particularly true because we will see that in the GenEvA algorithm, it is actually easier 
and more efficient to include leading-logarithmic resummation in matrix elements than not. 
However, we emphasize that the GenEvA algorithm can be used as phase-pace generator 
independently of any particular strategy for event generation, and similarly the GenEvA 
framework could work with an alternative phase space algorithm. 



2. OVERVIEW OF THE GENEVA ALGORITHM 

The idea of the GenEvA algorithm is straightforward: Generate events using a parton 
shower and reweight these events to any desired distribution in a way that avoids double- 
counting. In this section, we develop the master formula for the event weight and then 
summarize the various ingredients used in this master formula. The details of the algorithm 
are then given in Sees. [|, ||, and ||. Readers not interested in the technical details should 
get a basic understanding of the algorithmic issues involved from this section, and can then 
safely skip to Sec. [7] for a comparison to other event generators and the efficiency of the 
algorithm. 



A. The Master Formula 

Let $ denote a point in Lorentz-invariant phase space generalized to include all additional 
particle information, like flavor, helicity, and color, that fully characterize an event. We 
want to distribute events according to a general distribution cr($). That is, we want the 
contribution of an event with value $ to be da = er($) d$. Even though the value of cr($) 
is calculable for fixed values of $, it is usually very complicated to directly distribute events 
according to cr($). 

If we have a known normalized probability distribution for which it is easy to 

distribute according to, we can choose values of $ according to V(&) and assign the weight 

v ; P($) v ; 

to each event with value <3>, as in Fig. [I]. Since the probability of choosing a value $ is 
dV = V($) d$, the weighted set of events with values $ will be distributed according to 

da = w($)dV = ^^7>($) d$ = <r($) d$ , (2) 

as desired. 

Now imagine we have an easy way to distribute events in terms of a different variable 
X according to a known normalized distribution 'P(S), and there is a mapping $ = $(£) 
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FIG. 1: The simplest version of reweighting. Events are distributed according to a normalized prob- 
ability distribution V($>) and given weights w($) = a($>)/V($) such that the resulting weighted 
events are effectively distributed according to <r(3>). In GenEvA, corresponds to Lorentz-invariant 
phase space variables and cr($) is a differential cross section. 



P(£)J(£) 




FIG. 2: Reweighting using a different distribution variable E, where there is a one-to-one and 
onto map E — ► 5>(E). Events are distributed according to a normalized probability distribution 
"P(E), which together with the Jacobian J(E) from the variable transformation E — ► ^(S) defines 
an effective probability distribution V($>) = 7 ? (E)J(E). The weights can now be defined 

analogously to Fig. |l|. 



which is one-to-one and onto, as in Fig. Using "P(S) as sampling function, the probability 
for choosing a value of $ is now 



dV = V(E) dE = V(E) J(E) d$ = V{$) d$ . 
where E = S($) is the inverse of $(E) and 



1 



J(E) 



9(E) 



(3) 



(4) 



is the Jacobian in the transformation from E to $. Thus, we obtain the event weight 



P[E($)]J[E($)] 



(5) 



S 




FIG. 3: Reweighting when there are multiple values Ej that map to the same point $(Ej). If 
all Ej that map to a given $ are known, then there is again an effective probability distribution 
V($) = J2i "P(Ej)^(Ej) that can be used to define an appropriate weight w($>). In GenEvA, 
the parton shower provides a ^(E), and it is well-known that multiple parton shower histories 
Ej can correspond to the same phase space point 3>(Ej). For this reweighting strategy to be 
computationally feasible, there must be an efficient way to calculate the sum over all parton shower 
histories Ei^( s i)^( s i)- 



It is important that the mapping $(E) is onto, otherwise there would be values of <3? which 
are never covered by V(H). On the other hand, $(E) need not be one-to-one. Instead, for 
a given value of $ there could be a discrete number n($) of E values which all map to that 
same $ as in Fig. []. In this case, the probability of choosing a value $ is 

n($) 

dV = J2v(X i )J(X i )d®, (6) 

i=l 

where here and in the following Ej = Ej($) denotes the i-th point £ for which $(£j) = $ 
and i runs from 1 to n($). The weight as a function of $ is thus 

W$) = ^ (7) 

In GenEvA, the parton shower provides the sampling function V(H), where E describes 
a complete parton shower history, and $(E) is given by the shower's final state. It is 
well-known that multiple parton shower histories Ej can produce the same phase space 
point $(Ej). In practice, it may be computationally expensive to calculate the sum in the 
denominator of Eq. (0), as it requires reconstructing all possible parton shower histories Ej 
for a given point However, Eq. (|7|) is dictated by our choice of characterizing events by 
$, which forces us to treat the weight of an event as a function of $ as well. In particular, 
in Eq. (^) cannot depend on the actual Ej which was used to generate the event. 

To relax this condition we can simply characterize events by E and let w = w(E) be a 
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er($)d(Ei) 





cr($)a(S 2 ) 




FIG. 4: The reweighting approach used by GenEvA when the map <3?(S) is not one-to-one. Instead 
of taking the event weight to depend on only, the event weight w(E) is effectively a function of 
S. Using the overcounting factor the desired distribution <r($) is split up among the various 

Sj that map to the same including the Jacobian factor J(Ej). As long as ^i^K^O = 1, the 
resulting events will be distributed according to cr($) without the need to calculate 7 ? (Si)J(Sj) 
for every Sj. In practice, there is a small computational cost to make sure that a(£j) is properly 
normalized. 



function of E as in Fig. |4|. In this case, the total contribution to da for a given $ is 

da = J2 M^) dP, = £ ^ryP(^) dS = ^) dE . (8) 

At this point, we need to define what we mean by cr(E). From Eq. (Q) we have to require 

d<7 = 5^<7(E i )J(E < )d$ = «7($)d$. (9) 

i 

Note that simply taking <r(E) = c[$(E)]/J(E) does not work, because it would lead to an 
overcounting by a factor of n($) in Eqs. (pi) and (El). Instead, we need 



such that 



<r[$(E)]a(E) , v^^v,n 

^(S) = 1 l J( ^ } 1 ^ with & ^ = 1 



da = Y^ o-(Ei) J(S<) d$ = tr($) d(E;) d$ = tr($) d$ , 



(10) 
(11) 
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as required. 

Using Eq. (|l0|), the final weight as a function of S is 



Eq. (12) is the master formula underlying the GenEvA algorithm. Alternatively, if we wish to 
specify <r(£) directly, we use the notation o~i(&), where i labels the different shower histories 
mapping to $, and the event weight is 

It is clear that all the information in o"j($) could be absorbed into an ordinary cx($) times an 
overcounting factor However, note that the Jacobian </(£) still appears in Eq. (|13|). 

meaning that is still a differential function of $ and not E. 

B. Calculating the Event Weight 

The master formula in Eq. (|l2|) contains several ingredients. Here, we will give a brief 
overview of each of the components, while the details are discussed in the remainder of the 
paper. 



1. The Parton Shower Probability 'P(X) 

In the GenEvA algorithm, the parton shower plays the role of the sampling function 'P(S) 
for reweighting to any desired matrix element cr($). The crucial point is that the par- 
ton shower automatically contains the right singularity structure of QCD matrix elements, 
which should in principle allow a reweighting with reasonable efficiency However, for a 
parton shower to be used for the purposes of reweighting, it has to satisfy several additional 
requirements. 

First, we have to know the parton shower's probability distribution "P(S), which means 
we must be able to compute the exact probability with which the parton shower generates 



a given shower history E. The analytic parton shower algorithm constructed in Ref. |p2 
has this property and is the algorithm used here and detailed in Sec. [3|. The key feature 
of this shower is that energy-momentum is conserved by construction without any ad hoc 
momentum shuffling. 

Second, the parton shower has to cover all of phase space. Of course, almost any parton 
shower already covers all possible final state configurations, in that it will generate 2-body 
final states, 3-body final states, and so on, with the flavors of the final state particles 
representing all possible allowed configurations. What is important is that the shower also 
has no kinematic dead zones, to ensure that the map is indeed onto and all of phase 

space is covered as necessary. 

Finally, as already mentioned, the fact that the parton shower generates many final 
state particles seems at first glance like a liability, because matrix elements for ra-body 
configurations with arbitrary n are not readily available. In order to reweight the generated 
events, we need to be able to "backtrack" (not necessarily uniquely) a high-multiplicity final 
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state to a lower-multiplicity intermediate state, while still being able calculate analytically 
the probability with which this intermediate state was generated. We call this feature 



"truncation" , and it is possible because of the analytic properties of the algorithm of Ref. |52 
in conjunction with the parton shower's built-in notion of probability preserving ordering. 
This truncation is a concrete implementation of the phase space projection required by a 
phase space generator suitable for the GenEvA framework as explained in more detail in 
Ref. 0. 



2. The Jacobian J(S) 

An event produced by GenEvA's parton shower is characterized by a parton shower history 
E that consists of the set of 1 — > 2 splittings that build up the event. In contrast, Lorentz- 
invariant phase space $ characterizes an event by the set of three-momenta (symmetrized 
appropriately with respect to identical particles) of the on-shell final state particles. As 
needed for reweighting, there is a map from parton shower histories X to Lorentz-invariant 
phase space $, such that for every parton shower history £ there is a unique point $(£). A 
nontrivial Jacobian factor is required in changing variables from the natural variables 
used to describe the 1 — > 2 splittings in the shower history to those that describe Lorentz- 
invariant phase space. This Jacobian will be given in Sec. [|. 



3. The Overcounting Factor ce(S) 

It is well known that the parton shower does not map out phase space in a one-to-one 
way, because in general, multiple parton shower histories can lead to the same final state 
four- momenta. The overcounting factor takes care of this overcounting, and essentially 
determines how cr($) is split up among the different £«($) that map to the same point $. 
Since we are eventually only interested in o~(§), we can choose freely only subject to 
the constraint in Eq. (|lOl). The simplest possible choice would be = l/n($), where 

n($) counts the number of Ej mapping to the same $. This choice will in general lead 
to a poor statistical efficiency, since the weights will vary widely for the different £j 

contributing to the same point $. To ensure well-behaved statistics, we need to ensure that 
for a given phase space point $, all event weights are roughly the same, thus requiring that 
a(S) scales like J(E). We will discuss the choice of a(£) in the context of GenEvA in 
detail in Sec. |5|. 



4. The Matrix Element <r[$(£)] 

GenEvA can distribute events according to nearly any distribution cr($), and this quantity 
will depend on the underlying physics one is trying to describe. Crucially, physics consider- 
ations alone can determine the best choice for cx($), and the algorithmic details of GenEvA 
need not influence this choice. 

While o"(<£>) in Eq. ([12]) only depends on a set of on-shell four-momenta, in many cases 



it can be advantageous to use Eq. (13) and directly work in terms of o"i($). That is, we 
may want to assign different weights to events having different parton shower histories £j 
even though they contribute to the same point $(£«). Two reasons to use o"j($) are to 
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improve the efficiency of GenEvA and to resolve logarithmic ambiguities |[J]. In Sec. |6], we 
will discuss tree- level fixed-order QCD matrix elements using cr(<E>), as well as calculations 
which properly sum the leading- logarithmic dependence, using 0"i($). 



C. Reweighting Efficiency 



To assess the practical usefulness of the GenEvA algorithm, it is important to define what 
we mean by efficiency. In this paper, we will use statistical efficiency as the figure of merit 
in contrast to unweighting efficiency. Since unweighting efficiency depends on the value of 
the maximum weight generated, it is sensitive to how many events were generated in a given 
run. Here we argue that statistical efficiency is a much better measure of the strength of an 
algorithm. 

If TV is the total number of events, then in the limit of large TV the weights satisfy 3 

«)($) -> N /d$ . (14) 

n J 

To see this, consider the normalized distribution <t(<&) = cr($)/ fd$er($). Since both 
<5"($) and 7 ? ($) are normalized to unity, the normalized weights = a(§)/V(§) satisfy 
^ n w(<l ) ) — > TV, from which Eq. ( |14|) follows. 
Dividing Eq. ([14]) by iV we have 

_L = fd$a($), (15) 

n J 

where (•) denotes the sample average. From the standard error on averages, the statistical 
uncertainty on the Monte Carlo estimate of J d$er($) is thus given by 

^W-Wf , (16) 



N 

The numerator describes the width of the distribution of weights, which vanishes in the limit 
of all weights being equal. This is equivalent to the fact that for unweighted events, the 
uncertainty on the total number of events is always zero. Thus, the ratio 

_ M$)) 2 _ i [E,^)] 2 _iVe ff <1 (u) 

VCS (w(<S>) 2 ) TV £ n M$) 2 ~ N - ' ^ } 

provides a statistical measure of the efficiency of the reweighting. We can interpret N c q as 
the effective number of events after reweighting in the sense that the sample of reweighted 
events is statistically equivalent to an unweighted sample of TV e ff events. In the limit of 
all weights being equal ^ e ff — ► 1 and N e g — > N. Hence, this way of distributing events 



3 This definition of weights differs slightly from the one commonly used in Monte Carlo programs. Here, the 
weight of an event is independent of the number of events generated, whereas most Monte Carlo programs 
divide individual event weights by TV, the total number of events generated in a run. We will use the 
run-independent definition of weights throughout this paper. 



13 



by reweighting V(<&) to er($) is efficient as long as P($) is reasonably close to <x(<3>). For 
example, if we wish to distribute events according to a function cr($) which has large peaks, 
we have to choose a sampling function P(<3>) that has a similar peak structure. 

The definition of N eS differs from that of N unw , which is the number of events that would 
remain after unweighting the sample, 

6ff " £„^) 2 mm ~ maxK$)] ■ (18) 

The unweighting efficiency r] unw = N unw /N is highly sensitive to the tail of the weight 
distribution, such that even an exponentially small number of high-weight events can cause 
a dramatic decrease in r] unw . The only way to deal with this is to throw away events with too 
large weights, which of course introduces a systematic error which is very hard to quantify. 
In contrast, the statistical efficiency N e s/N correctly reflects the full statistical power of the 
event sample, and being based on sample averages, is unaffected by individual weights. 

The main disadvantage of considering the statistical efficiency is that since the statistical 
uncertainties are determined by iV eff and not by N, it is necessary to store and process a 
larger total number of events than in a statistically equivalent unweighted sample. However, 
if this becomes a serious issue one has always the option to partially unweight the sample 
by unweighting just those events that have small weights. That is, we can define a threshold 
weight Wo, and only unweight those events with weights smaller than wq. The expected 
statistical efficiency of the partially unweighted sample is given by 

Vcs{w ) = — ^ ^ 7-- • (19) 



Note that r/ e fr(0) = i] c g and rj e s(w max ) = rj unw . In Sec. [70. we will report the efficiency of 



GenEvA not only in terms of ^ e fj, but also in terms of the time needed to achieve a partially 
unweighted sample with rj e s{w ) = 0.9. 



3. THE ANALYTIC PARTON SHOWER 

As discussed in Sec. |2 B 1| , a parton shower to be used as phase space generator has to 
fulfill three requirements: 1) It has to have an analytically calculable probability "P(S); 2) 
it has to cover all of phase space; and 3) we have to be able to analytically truncate the 
shower, which means we must be able to truncate any high-multiplicity final state produced 
by the shower to an ra-body final state for any given n, while retaining analytic control over 
the truncation probability. In the following we will describe in detail how these requirements 
are met by GenEvA's parton shower. 



A. Overview 



The parton shower used by the GenEvA algorithm to generate phase space is based on 
the analytic parton shower algorithm constructed in Ref. [52]. As we will only be discussing 
e + e~ — > n jets, we will only discuss the final state parton shower. We comment on the 
extension to hadronic collisions and initial state showers in Sec. B A. 
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FIG. 5: The structure of the parton shower. A solid line indicates an already processed particle, 
while a dashed line indicates a particle that has not yet been processed. A dot at the end of 
a solid line indicates a particle that did not branch above the shower cutoff, i.e. a final state 
particle. Starting from an existing branch M — > LR, where a mother particle M has branched 
into daughter particles L and R, GenEvA's parton shower generates the double branch M — > LR — > 
(LL/LR)(RL/RR), where it is understood that there is also a finite probability for L or R to 
remain unbranched. Understanding how the L and R branches are coupled is essential for being 
able to exactly compute the double-branch probability and to truncate the parton shower. 



Since the main purpose of GenEvA's parton shower is to generate n-body phase space 
efficiently, the current implementation is only accurate to leading-logarithmic accuracy, and 
neglects several important physics effects a realistic parton shower should take into account, 
such as the running of a s and color coherence |5(J |73|]. In addition, the shower we consider 
here treats all final state particles as massless. While it is straightforward to construct an 
algorithm that does incorporate all of these physics issues and still maintains analytic control 



(see the discussion in Ref. ||52||), for the purposes of this paper the simple analytic parton 
shower discussed here will suffice. 

An event produced by GenEvA's parton shower is characterized by the kinematics of the 
initial hard-scattering process (typically a 2 — > 2 interaction), followed by a set of 1 — > 2 
particle splittings. We mostly follow the notation of Ref. . We call a given 1 — > 2 splitting 
a branch, with the original particle called the mother particle, and the two final particles 
called the daughter particles. We label the mother particle by M, and the daughter particles 
by L and R. The daughters of these daughters (the granddaughters) are labeled by LL, LR, 
RL, and RR. 

Starting from an existing branch M — *> LR, the shower generates a double branch M — > 
LR — > (LL/LR)(RL/RR). This is illustrated in Fig. [|, where here and in the following 
figures a solid line indicates an already processed particle, while a dashed line indicates a 
particle that has not yet been processed. A dot at the end of a solid line indicates a particle 
that did not branch above the shower cutoff, i.e. a final state particle. The double-branch 
probability, i.e. the probability to obtain a specific 1 — > 2 — > 4 double branch from a given 
1 — > 2 single branch, is the fundamental object in GenEvA's parton shower, and allows it to 
have a well-defined notion of global evolution. 



15 



L 




FIG. 6: The kinematic variables describing a single splitting in the parton shower. The splitting 
functions can be chosen to be functions of the mother's virtuality tM and the angles cos 9m and 
4>m, where cos 9m is the angle between the L daughter and the boost axis to the mother's rest 
frame, and 4>m is the azimuth around the boost axis. Note that the L and R four-momenta cannot 
be determined until ti and tn are known, that is, only after the daughters themselves have been 
processed. 



B. Choice of Kinematic Variables 

In general, the splitting functions describing the 1 — ► 2 parton splittings depend on 
three kinematic variables, whose precise definitions depend on the details of the parton 
shower. One of these variables is the evolution variable of the parton shower, which is always 
decreasing as the shower progresses, while the other two describe the remaining kinematics 
of the splitting. There are several choices possible for the evolution variables, such as 
virtuality |5J| the energy- weighted emission angle 0, [7i| , transverse momentum |TUL [751, 
, or more general variables []33 



In GenEvA we use the analytic parton shower algorithm denned in Ref. |5JJ, in which the 
evolution variable is the virtuality tM of the mother, and the remaining two variables are 
chosen to be the angle cos 9m in the mother's rest frame between the L daughter and the 
boost axis, and the azimuthal angle (j>M around this boost £1X1 S ; cLS shown in Fig. || Thus, 
we can write the differential kinematics of the parton shower in terms of 

dS 2 (M -> LR) = dt M dcos0 M d0 M , (20) 

with the phase space limits 

-1<COS0 M <1, O<0Af<27T, 0<VtL + VtR<VtM- (21) 

The primary reason to use the angles as the basic variables is that their phase space limits 
are trivial and completely independent of any other kinematic quantity, which also makes it 
trivial to satisfy energy-momentum conservation for each branch. 

A complete shower history with n final state partons is described by 

n-2 

dS n = d$f n dS 2 (M, -> LM , (22) 

i=l 

where is the two-body phase space of the initial hard scattering and dS 2 (Mj — > LiRi) 
are the subsequent branches from the shower. We will discuss the precise phase space 
measure and its relation to ordinary n-body phase space in more detail in Sec. §. 
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Instead of the angle cos 9 m, many parton showers use the energy splitting zm = El/Em, 
describing how the energy of the mother is shared between the two daughters. Indeed, the 
splitting functions themselves are usually written in terms of z. The relation between z and 
cos 9 is given by 



cos9 M = -t— 



Z M = ~ 



1_ %Zm 

@M 
1 

2 



[l + t L /t M -t R /t 



Mi 



A(iM) ti, tit) 



1 + j— — + /3m cos 9 m A(t A /; t L , t R ) 
im £m 



(23) 



where tM is the virtuality of the mother, ti,R are the virtualities of the daughters, and we 
have defined 



P. 



t 



M 



M 



E 2 

M 



A(*m; t L , in) 



- — a/ (tM — t~L 

tM 



t R ) 2 ~4:t L t R . 



(24) 



with Em being the energy of the mother. 

In GenEvA's shower, cos 9 is the fundamental splitting variable, while z is a derived quan- 
tity. However, since the usual splitting functions are functions of z, we need to translate cos 9 
into z. At the time the mother is branched, both daughters are still massless, so Eq. ( PB"| ) 
gives the relation 

z(cos 9m, E m ) = - + ~Pm(E m ) cos 9 m , (25) 

where we made the dependence on Em explicit. The problem is that the final value for the 
mother's energy Em cannot be determined until after the mother itself and its sister have 
been processed and their final invariant masses are known. This is equivalent to the fact that 
in ordinary parton showers based on 1 — > 2 splittings, the value of zm has to be "re-shuffled" 
at a later stage in the algorithm in order to satisfy momentum conservation. For a more 
detailed discussion see Ref. 



52 



In contrast, once a value for cos 9m is determined, it never 
has to be changed afterwards, which crucially relies on the fact that its phase space limits in 
Eq. (|2~l"D are independent from the rest of the shower. Nevertheless, we still need to pick a 
value of Em in Eq. ([25]) , which we denote as Efl rt . From Eq. (|25|) we see that the choice of 
Efl rt effectively determines the range of z allowed in the splitting functions. Therefore, to 
avoid artificially large subleading-logarithmic effects, one should pick an E s ^ rt of the order 
of Em- In our case, the precise choice of Ef^^ is dictated by the requirement to have a 
simple analytic truncation, and will be discussed in Sec. |3 E| . 

We are now ready to obtain the required splitting functions. Taking into account the 
appropriate Jacobian factors, we have 



f M ^LR(tM, cos 9 



m, <pm; E m 



start \ 



? % ^ ^~ Pm^lrHcos 9 m , Etr')} 



2ty 2 2tt tM 

where Pm^lr( z ) are the well-known Altarelli-Parisi 0] splitting functions 

1 + z 2 



(26) 



Pq^qg(z) — C F 



1 - z 
1-z 



+ 



z 



1 - z 

P g ^(z) = T R [z 2 + (1 - zf] 



+ z(l-z) 



(27) 
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FIG. 7: The single-branch probability. Starting at the scale £ s tart> the mother M can either branch 
at the scale tM > tent or remain unbranched down to the scale t cu t, which is denoted by a dot at 
the end of a line. The branching probability V^^Mi vm) an d the no-branching probability Vff 
also depend on i s tart, *cut, and Eff vt . The variable vm is a shorthand for the splitting angles and 
the daughters' flavors. 



In the following, to simplify the notation, we will mostly suppress the dependence of the 
splitting functions on Ef;^, but one should keep in mind that wherever the splitting func- 
tions appear there is also an implicit dependence on E s M art . 



C. The Single-Branch Probability 

The basic ingredient in the parton shower is the single-branch probability, which deter- 
mines if and how a given mother particle M splits into two daughters L and R, as illustrated 
in Fig. [7]. The probability for a single branch to occur with given values of the splitting vari- 
ables is a combination of the splitting function and a Sudakov factor (for a review of parton 
showers see for example Ref. ||): 

'Pm-^LR^M, COS 9 m, 4>M] istart, ^cut) = /m^LK^M, COS 9 m-, 4>m) A-M^startj thl) 

X 9(t sta ,rt > tf^) 9{t M > tout) ■ (28) 



The Sudakov factor is defined as 

-t 2 rl r2n 



A M (*2,*i) = exp 



pt2 l>i /*Z1T 

/ dt dcos# / def) fhi^xjt, cos 9 
Jti J-i Jo x 



(29) 



where X denotes all possible flavor structures allowed for the splitting of the mother M. 
The Sudakov factor can be interpreted as a no-branching probability, since the probability 
of obtaining a value t M < t s tart is given by the probability to not branch anywhere between 
t m and t s tart times the probability to branch at tM, which is given by the splitting function 
evaluated at tM- Thus the probability of having no branching anywhere above t cut , where 
t cut is the cutoff scale of the parton shower, is determined by the no-branching probability 
directly 

(^startj ^cut) = AM^start; ^cut) • (30) 

For simplicity, we will frequently use the shorthand notation 

Vm^Mi VM] ^start, ^cut) = ^M^LR^M, COS 9 Mi 4>M] ^start, ^cut) , (31) 
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where vm denotes all variables other than the virtuality required to describe a single splitting. 
Thus, vm contains the angular variables cos 9 m and <Pm, as well as the flavor information 
M — > LR. We will also use the notation 



1 p2tt 



1 JO 



dv M = / dcos0 / d0 ■ ( 32 ) 



X 



Of course, the probability for a particle to do anything (branch at any tM > t cut or not 
branch above t cut ) is equal to unity 

rtstmt r 

/ d^M / dt?M Vm^m, Vm] t star t,t cut ) + V M (tstart, tcut) = 1- (33) 

•'tcut J 

To see this explicitly, note that the Sudakov factor satisfies 

-A M (t 2 ,ti) = A M (t 2 ,ti) / dcos6» / d0 fM->x(h, cos 0, 0) 
f l 7-1 Jo x 

= J dv M V^{t u v M ] t 2 , tcut) , (34) 



d 

d^ 



where we are assuming that ti > t cut . This implies 

dtjif J dVM Tm^MjVm'i tstartjtcut) = AM^start) h) ~ A^f (tstart; ti) , (35) 

from which Eq. (|33"D follows. 

D. The Double-Branch Probability 

GenEvA's parton shower algorithm always acts on an existing M — > LR branch. It takes 
both the L and R daughters of the branch and considers their branching into two sets of 
granddaughters LL/LR and RL/RR. The shower algorithm is illustrated in Fig. |8], and 
works as follows: 

1. Pick a random branch with two unprocessed daughters L and R and a mother particle 
M with invariant mass £m- 

2. Choose Ef 3 ^, set t start = t M , and determine values of {tx,fi} and Vl} indepen- 
dently according to the single-branch probability V\ x {Ui v i] tM,t cut ), where i = L,R. 
If ti < t cu t, leave that particle unbranched. 

3. Keep the branch of the daughter with the larger t, call this the max-daughter, and 
define t max = max(t£,, Discard the branch of the min-daughter. 



4. Determine new values {t m i n , f mm } for the min-daughter according to 

^min(tmin,^min; t*,t cut ), with t start Set to 

t„ = min [t max , (y/tM~ - Vtmax) 2 ] • (36) 
If tmin < tcut, leave that particle unbranched. 
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FIG. 8: The GenEvA shower algorithm. Starting at the scale tM, trial L and R branches are 
generated, but only the one with the larger of and tR is kept. Then, starting from the new scale 
i* = min[t max , (\/Tm — vt) 2 ] , the other branch is generated. This coupling between the L and 
R branches is essential to obtain an analytic truncation, and means that the GenEvA shower has a 
notion of global evolution. 



5. Continue until all branches are processed. 

Since Eq. ( pq) contains no restrictions on the angles cos# and 0, their phase space is 
trivially covered. The range of tL,R allowed by the algorithm is 

£cut < tl,R < tM , Vtl + VtR < \ft~M~ , (37) 

where the second condition is enforced by t* in step 4 of the algorithm, as required by 
Eq. ([21]) . Hence, apart from the infrared shower cutoff t cut , the algorithm covers all of phase 
space without leaving any dead zones. In particular, a parton shower history with a cutoff 
value t cut and n particles in the final state covers the entire n-body phase space with the 
additional restriction that the minimum virtuality of any two particles is greater than t cut . A 
As there are QCD singularities associated with two partons getting close in virtuality, the 
t cut restriction is actually helpful to regulate infrared divergences in tree-level diagrams. 

The above algorithm results in an analytically calculable parton shower. It is clear that 
the algorithm couples the splitting of the L and the R daughters together, resulting in a 
double-branch probability 

VM{tL,VL] t R} VR] tM,t C ut) =VM^LR^(LL/LR)/(RL/RR)(tL,COs8 L ,(f) L ; t_R, COS 9r, (j)R] t M ,t cut ), 

(38) 

where t^R are the final virtualities selected for the daughters, and vl,r denote all other 
variables describing the daughters' splittings, including the flavor and color structure for the 
LL, LR, RL, and RR granddaughters. 

There are four possible cases for the double-branch probability as shown in Fig. ^ 1) 
Both the left and right daughters branch; 2) the left daughter branches and the right one 
does not; 3) the right daughter branches and the left one does not; and 4) neither the left nor 



4 Two partons from different branches can also accidentally get closer than i cut , so there are additional 
regions of n-body phase space that are covered in practice. We will discuss this issue more in Sec. |6 A|. 
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FIG. 9: The double-branch probability. The GenEvA algorithm in Fig. || yields four different 
branching/no-branching combinations. They are built out of products of single-branch probabili- 
ties, with nontrivial interdependence between ti and t R . 



the right daughter branches. Following the shower algorithm, the double-branch probability 
for case 1) is given by 

^M br (^> V L] t R ,V Rl ^M,^cut) 



0(t L > t R ) V b L r (t L ,v L ; t M ,t cut ) x A R {t M ,t L )P^(t R ,v R ; t*,t 



cut; 



+ 0(t R > t L ) [V R T (t R , v R ; t M , t cut ) x A L (t M , t R ) Vl T (t L , v L ; t*, t cu t)J , (39) 

where the extra factors of Al ±r arise from discarding the first attempt to determine the 
smaller t value. For the other three cases we have 



(tL) V Li tM,t cnt ) — V L T (t L} VL, tjv/,t cut ) x A R (tM,t L ) A R (t*,t, 



cut) i 



T^M T (tR,V R , tM,t cut ) — V R (t R ,V R ] tM,t cu t) X Al^Mi tit) AL(t*, tcut) 



V M 11 (tfcfj tcut) — A L (t M ,t cut ) X A R (tM,t, 



cut ) 1 



(40) 



where it is understood that if t* < t cut , then V™{ti,Vi\ t*,t 



CUt ! 



and Aj(t*, t. 



cut ; 



We can formally combine all four cases by writing 

Vm^LiVl, t R ,V R ] tMjtcut) 

= *Pm T (tL, Vl] t R , V R , t Ml tcut) 0(t L > t cut ) 9(t R > t cut ) 

+ ^M nb (tL, v L ; t M , t cut ) 6(t L > t cut ) 6(t R ) 5(v R ) 

+ ^M hT {t R) v R ; t M , tcut) 0(t R > t cut ) 8(t L ) 8(v L ) 

+ ^ f bnb (t M , tcut) 8(t L ) 5(v L ) 8(t R ) 8(v R ) . (41) 

The 5 functions are just formal crutches to be able to define values for tL,R and vl jR when no 
branching actually occurred. When a daughter does not branch, it is put on shell through 
the 8(tL,R) functions (for massless particles as we are considering here). The 9 functions are 
not strictly necessary, as they are already contained in the single-branch probabilities, but 
to be explicit we write them out again. 

As a check of the algorithm, it is straightforward to verify that the probability for anything 
to happen is equal to unity: 



/tu r ftM r 

dt L J dv L J dt R J dv R V M (t L ,v L ; t R ,v R ; t M ,t cut ) 



1 . 



(42) 
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FIG. 10: An example parton shower history with 5 particles in the final state. The quantity 
d<7Af^£ij/<T t otai is the probability to obtain the initial M — > LR branch from a 2 — > 2 hard 
scattering as a function of $2- The total probability for obtaining the full 5-particle final state is 
given in Eq. fl44|). 

To see this, integrate over min(tL,t R ) using Eq. (|35f) , and then examine the quantity 

^[A L (t M ,t)A R {t M ,t)] , (43) 

to obtain a simple expression for the integral over max(tL,t R ). 

It is now straightforward to calculate the probability for generating an event with a given 
hard 2 — > 2 scattering and a history of branches with n particles in the final state. It 
is simply given by the probability to obtain the hard scattering multiplied by the double- 
branch probabilities for each of the following branches. As an example, consider the parton 
shower history shown in Fig. pi]. The differential probability to obtain this event is given by 

dP(S Fig .|i^D = d ° M ^ LR x p^ hr (t L ,v L ; t R ,v R ; t M ,t cut ) dt L dv L dt R dv R 

x Vl Tnh (t LL ,v LL ; t L ,t cut )dt LL dv LL x P£ bnb (t fi ,t cut ) 

x P£ b L nb 0Wcut) , (44) 

where do m->lrI 'ctotai is the probability to obtain the initial M — > LR branch from a 2-body 
matrix element with L and R in the final state. Hence, GenEvA's parton shower algorithm 
has indeed an analytic formula for the probability V(H) to obtain a given parton shower 
history. 



E. Truncation 

The parton shower described above will produce final states with in principle arbitrarily 
many final state particles. The goal of GenEvA is to reweight such an obtained final state 
to a known distribution. In practice, the exact distributions will only be available for 
a limited number of final state particles. Furthermore, for the purposes of the GenEvA 
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framework [l| , one might want to choose the scale that separates partonic calculations from 
phenomenological models to be different from the value t cut used by the parton shower. 

Consider the case where we only have partonic calculations available with up to n max 
partons in the final state, and we want the minimum virtuality described by the partonic 
calculation to be t ma tch > tcut- The parton shower algorithm will in general produce events 
with more than n max partons in the final state and produce splittings with t < t ma tch- By 
assumption, these events cannot be reweighted to the partonic calculation. One way to 
deal with such events is to simply reject them, and only reweight events that satisfy the 
criteria required by the partonic calculation. However, this would reject a large number of 
events and would obviously result in large inefficiencies. Furthermore, because the parton 
shower preserves probability, unless the user calculated the rejection probability, vetoing 
events would make it impossible to obtain total cross section information from the phase 
space generator. 

Instead of rejecting events with too many final state partons or too low a splitting scale, 
we can truncate these events to a final state with at most n max particles, in such a way that 
each branch has virtuality above t ma tch- This is achieved by a simple truncation algorithm 
that truncates events to a given {t ma tch, n miB }: 

1. Order the branches by the value of the evolution variable t. 

2. Remove any branches with t < t ma tch- That is, remove the daughter particles from 
any mother whose virtuality is t < t ma tch, leaving the mother unbranched. 

3. If the final state contains more than n max particles remove additional branches, starting 
with the smallest value of t, until there are n max final state particles. 

4. Recalculate the kinematics of all remaining final state particles. 

5. For events with n max final state particles, redefine t ma tch to be the smallest value of 
the virtualities t of the remaining branches. 

This algorithm clearly results in an event sample for which all branches have t > t ma t c h 
and there are no more than n max final states. The last step guarantees that if the shower 
restarts from the event-defined t ma tch scale, then the resulting distribution of events will 
be identical to those if truncation were not applied. Note that this algorithm preserves 
probability because no events are ever rejected. 

To be able to reweight the resulting events after truncation, we of course need the proba- 
bility with which the truncated event was generated. As expected from the notion of global 
evolution, the probability for finding a truncated event should be exactly equal to the prob- 
ability of generating that event with the shower cutoff t cu t set equal to the event-defined 

tmatch; 2.6. 

P^ UnC (tL, Vl, tji,VR] t M ,t cut ] tmatch) = Vm^L^L\ tR,VR] t M , t ma tch) . (45) 

To see this, we write out V^ nc in analogy with Eq. ( pTTp as 

V^ nC (tL,Vi; tft, Vr] tM,t C nt] t ma tch) 

_ -ptruncbrbr^^ f R ^y R ^ t M) t cut ; t ma tch) 0(t L > t m atch) 0(t R > t m atch) 
+ ^M UnCbrnb (t,L, Vl] t M , t cut ; tmatch) 9(tl > t m atch) ^(ti?) S(Vr) 
+ ^M UnCnbbr (^, Vr] tM, tcut! t ma tch) 0(t R > t ma tch) S^l) 5(v L ) 

+ K uncnhnh (hi, t cut ; Wh) S(t L ) 5{v L ) 5(t R ) 5{v R ) , (46) 
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FIG. 11: To get a truncated event where both daughters are branched, the original event has to 
have t LjR > imatch- 

AFTER TRUNCATION BEFORE TRUNCATION 




^ ^match ^match ^R ^ ^match ^match 

FIG. 12: To get a truncated event where only the left daughter is branched, the original event has 
to have > i ma tch and tR < t matc h- Since, formally, an unbranched particle has t < t cut , there are 
two different ways this can happen. For the truncated event, the shower cutoff is effectively raised 
from t cu t to t m atch- The same holds if the roles of L and R are interchanged. 

and show that 

M ' ' i ^Mi ^cut! ^match) = ^ M {' ' ' j ^Mj ^match) • (47) 

There are three cases to consider. First, take a double branch for which truncation keeps 
both the left and the right branches. This occurs if t L and tR are initially greater than t matcn , 
as in Fig. |TT[ Thus, the probability of getting given and tR values after truncation is 

ptruncbrbr^^ f R ^ VR ^ t Mj t cut ] t matc h) = '(^Ij ^Li ^Ri^r] ^Mj^cut) 

= V^ hr (t L , v L ; t R , v R - t M , t match ) , (48) 

where we used that the single-branch probability in Eq. (|28| ) only depends on the shower 
cutoff through an overall 9 function. Thus, the probability to get this double branch after 
truncation is equal to the probability to get this double branch from running the shower 

With t cut = t ma tch- 

Next, take the case where after truncation one daughter is branched and one is un- 
branched. This occurs for example when tR < t matc h < tz, which happens if the left daugh- 
ter branched above t m atch, while the right daughter did not branch at all or branched below 
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t R < imatch 

FIG. 13: To get a truncated event where neither daughter branches, the original event must have 
tL,R < ^match- Since, formally, an unbranched particle has t < t cu t, there are now four different 
ways for this to happen. For the truncated event, the shower cutoff is effectively raised from i cut 

to imatch' 



imatch, as in Fig. [L^. In this case, the truncation probability is 

itrunc br nb 

'""match 



-ptruncDrriD (f L ^y L ^ t M ,t cut ; t matc h; 



feu 
br/ 



dt R I dvRV h ^ bl (t L ,v L ; t R ,v R - t M , t cut ) + V$ nb (t L , v L ; t M ,t cut ) 



br nb / 



dt R / dv R V R (t R ,v R ; t„,t cut ) +A R (U,t c 



= V L \t L ,v L ; t M ,t cut ) x A R (t M ,t L ) 

= Vi r (tL,VL', tM,£match) X Afl(t&f, t£) A R (t*, t ma tch) 

= V% nb (t L ,v L ; t M ,t match ), (49) 

where we used the definitions in Eqs. fl39|) and p0|), and in the third step we used Eq. (|35| ) and 
again that V br only depends on the shower cutoff through a 9 function. The crucial feature 
of GenEvA's shower algorithm that makes this calculation possible is that the double-branch 
probability can be integrated analytically over the smaller virtuality t R without knowing the 
value of the larger virtuality t L . An analogous calculation holds for < t match < t R . 

Finally, take the case where after truncation, both daughters are unbranched. This occurs 



if both daughters either did not branch or branched below t ma tch ; as shown in Fig. 13. In 
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this case, the truncation probability is 



T^trunc nbnb /i . . ± s 
I M y^Mi tcut; t-matchj 



tmatch P /"^match P 

dt L / dv L / dt R / dv R VM hl {t L ,v L ; t R ,v R ; t M ,t cut ) 

tcut J J tcut ^ 

/"^match /* 

+ / dt L dv L V^ nb (t L ,v L ;t M ,t cnt ) 

''^mat ch 



b CUt } 



match P 

+ dt R dv R V^ hr (t R} v R ; t M , t cnt ) + V$ nh (t M , t c 

J tcut J 

_ pnbnb^^ t match ) , (50) 

where we used Eq. ( f42| ) and the fact that the truncation algorithm preserves probability. 
One can also do the ti and t R integrals explicitly as in Eq. ((49|) to verify this result. The 
above calculations show that Eq. ( [25]) indeed holds: The probability to get a truncated event 
is exactly the same as the probability for it to be produced by the shower with t cut = t ma tch, 
consistent with the notion of a global evolution. 

There is one subtlety we have completely ignored so far. In step 2 of the shower algorithm 
outlined at the beginning of Sec. [3 D| , we need to choose values for Ef; 3 ^ to be used in the 
daughters' splitting functions. To not violate Eq. (|45|), we need a method to determine 
values for Ef 3 ^ which gives the same values before and after truncation. In other words, 
we need an expression for E^ff that is invariant under truncation. 

Certainly the easiest choice would be to take Ef^ = Eqm- However, as noted in Sec. |3B| , 
we would like Ef; 3 ^ to be of the order of the daughters' actual energies. In principle, there are 
many possible choices. In analogy with t start = t M , we might attempt to take Ef 3 ^ = E M , 
where Em is the mother's energy. However, this does not work, because Em is not invariant 
under truncation. 5 Instead, we take Ef^ 3 ^ = E M 3 *, where E M 3X is the maximum possible 
energy of the mother and is determined recursively as follows: Given the Effi* used to branch 
the daughters, we define the daughters' E^ 3 ^, to be used in the subsequent branching of the 
granddaughters, by 



rrimax rpinax 
h L,R = 



l + ^+^ncos^f 1 

tM V 



M v t M 



(51) 



This corresponds to taking the minimum and maximum values of Zm in Eq. ([23|). For 
the initial hard scattering we simply take E M 3X = -Ecm- This definition of £J max is indeed 
invariant under truncation, since it only depends on E C m and the values of cos 9 and t of a 
particle's direct ancestors, i. e. only on those branches that cannot get truncated unless the 
particle itself is truncated first. 

With truncation, we now have two different ways of using the GenEvA algorithm in the 
context of a complete event generation framework, as illustrated in Fig. pj]. One method 
is to run the GenEvA algorithm to t cut ^> Aj R and truncate it to {t CVLt , ^max}, which means 
that {n < n max )-parton events are unaffected, while (n > n max )-parton events get truncated 
t° ?Vax-parton phase space with a variable, event-specific t matc h scale. The events can then 



5 As shown by Eq. (^) , the energies of two sisters depend on each others virtualities as well as their mother's 
energy and cos 6. Thus, truncating the aunt of a particle changes its mother's energy. 
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FIG. 14: Two different ways of using the GenEvA algorithm in a complete event generation frame- 
work. In the left panel, the GenEvA algorithm is only used above a scale fi = \/t cut 3> Atr, and 
traditional showering and hadronization is applied starting at fi. In the right panel, GenEvA is used 
to distribute events all the way down to v^cut = Atr to interface directly with a hadronization 
scheme. The first case allows the use of already existing parton showers that have been tuned to 
data, while the second case allows the scales associated with the partonic calculations to vary even 
after hadronization and detector simulation. In the companion paper ||], we discuss how both 
methods can give the same leading-logarithmic Sudakov improvements. 

passed to an external showering/hadronization program to start a new shower at t cut for 
n < Wmax and the variable t ma tch scale for n = n max . This uses GenEvA more-or-less as a 
traditional phase space generator. Alternatively, we can run the GenEvA algorithm all the 
way down to t cnt = Af R , interface with a hadronization routine, and then after the fact 
truncate the shower history back up to {t ma tch 3> Aj R , n max } to obtain a point in n-body 
phase space which now already has a complete shower attached to it. 

Because truncation preserves probability, both of these methods will give statistically 
identical results, however, the second method allows multiple different truncation schemes 
on the same already hadronized (and even detector-simulated) event. Of course, in prac- 
tice an external showering program can be tuned to data, whereas the GenEvA shower only 
guarantees leading-logarithmic accuracy, but the possibility to improve perturbative theo- 
retical distributions after detector simulation is an intriguing possibility, and for this reason 
it might be useful to develop an analytically calculable but realistic parton shower. 



4. THE JACOBIAN FACTOR 

In the previous section, we saw that the analytic parton shower covers all of phase space. 
However, the parton shower distributes events in terms of the shower variables {t j cos 6, </>} as 
in Eq. (|20|) , and there is a Jacobian factor from converting these to the variables describing 
Lorentz-invariant phase space. 

The complete parton shower space dS n for a shower history with n final state particles 
has two pieces. First, there are the variables of the initial hard scattering that yields 

two on-shell (and for us massless) final states A and B. Second, there is the product of 
the individual dS 2 splitting variables that describe the subsequent 1 —>■ 2 parton shower 
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branches. This gives 



n-2 



d£ n = d^focMji^Pfl) U^M, - L^) , (52) 



i=i 



where are the four-momenta of the particles A, B before running the parton shower, 
and i runs over all parton shower branches. We define that the way the parton shower takes 
the momenta of the A and B particles off shell is by leaving the angles df2 2 unchanged. 
Therefore, the phase space for the initial on-shell final state particles is related to that of 
the ultimate off-shell, showered particles via 

&®T{Pcm,Pa ,P { b) = T7T \ . n d$ 2 (pcM;PA,Pg) , (53) 

a^icm; oa, OB) 



where the ordinary phase space factor A is given in Eq. 

To convert the natural shower variables dS n to Lorentz- invariant n-body phase space 
d$ n , we can use the well-known method to recursively decompose n-body phase space in 
terms of multiple copies of 2-body phase space: 

d$ n (p CM ;Pl, ...,Pn) = (27r) 4 5 4 (pcM - X>J P ff (P? - ^ fl (P*) 

dQ 2 

= d$ n _i(p C M;Pl, • • • ,Pn-2, Q) -7T—d®2(Q;Pn-l,Pn) 

111 

= d$ n _i(p C M;Pi, • • • ,Pn-2,Q) -jt; — — T dt Q dcos 6 Q > ( 54 ) 

■J \pQ] t n _l, l n ) 

where the 2-body Jacobian factor is given by 

64tt 3 

J{t M ;t L ,t R ) - — — - — — -. (55) 

A(£Af, CLj or) 

The variables £q = Q 2 , cosOq, and 0q in the last line of Eq. ([jj) are exactly the same vari- 
ables used in d£ 2 in Eq. (^0| ). Their integration ranges are such that for any d$2(pM',PL,PR) 
we have 

-l<COS0Af<l, O<0m<27T, < V / tL + V 7 ^ < V 7 ^ , (56) 

identical to Eq. (|21|). Lorentz-invariant phase space can thus be written suggestively as 



n-2 

d$ n (PcM;Pi, • • • ,p n ) = ^ 2 (Pcm;Pa,Pb) TT -777 - — — — d£ 2 (Af; -> L;#i) , (57) 

~ = "7 'J{o Mv t Lv t Ri ) 

where A and .B can be thought of as the particles originating from the hard interaction and 
tMi, and denote the virtualities of the mother and daughter particles in the i-th phase 
space recursion. 

Though both the parton shower and Lorentz-invariant phase space can be built re- 
cursively, there is one important difference between the parton shower dS n and Lorentz- 
invariant phase space d$ n regarding identical particles. A parton shower distinguishes all 
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final state particles by a series of "L" and u R n labels, regardless of whether those parti- 
cles have the same quantum numbers. In contrast, Lorentz-invariant phase space does not 
distinguish between identical particles. Thus, an additional factor has to be included in 
the Jacobian to account for this difference. For a parton shower history with n final state 
particles that contains ni, n 2 , . . . , n k sets of identical particles, with n\ + n 2 + ■ ■ ■ + n k = n, 
we have to divide the Jacobian by an extra factor of ni! n 2 \ • • -n^-- 

Putting all this information together, the total Jacobian relating n-body parton shower 
space to Lorentz-invariant n-body phase space is 



J(E r 



1 



n-2 



A(*cm;*a,*b) ni\n 2 \ ■■■n k \ 



| \ J it Mi ? thi , tRi ) 



(58) 



i=i 



This factor is a crucial piece of the master formula Eq. (IT 



5. THE OVERCOUNTING FACTOR 
A. Overview 

Different parton shower histories can give rise to the same set of final state four-momenta, 
so the parton shower covers phase space multiple times. Therefore, if we wish to reweight 
a parton shower event to a matrix element, which only depends on the phase space $ of 
the external particles but knows nothing about the underlying parton shower history £, we 
have to account for the fact that the mapping $(£) from parton shower space £ to phase 
space $ is not one-to-one. This is achieved by the overcounting factor d(£) in Eq. dl2|). 

Since there is only a fixed number of parton shower trees that can yield a given set of 
external particles, there is also only a fixed number n(<3>) of points Ej($) that map to a 
given point $. Certainly the easiest choice for d(E) is thus 

which only requires the knowledge of n($) for $ = $(E), but not the precise form of each 
E,;($). Furthermore, we have defined our parton shower without any dead zones, such that 
for a given set of n final state particles, each compatible parton shower tree fully covers 
the corresponding n-body phase space. Hence, n($) only depends on the particle types in 
$ and thus equals the number of compatible parton shower trees, a quantity that could 
be calculated in advance. On the other hand, the reweighting efficiency of this choice will 
usually be poor, because if two points Ei i2 (<1>) have V(Ei) ^> 'P(S 2 ) then Eq. fljjD yields 
weights u>(Ei) w(Y; 2 ) that are very different. 

The optimal choice in terms of reweighting efficiency is to have the same weight w[Ej($)] 
for all i and a given $. This is achieved by 

a(E) = E,P(£,)W (60) 



which just reduces Eq. (|T2|) back to Eq. (0). Unfortunately, calculating the sum in Eq. (j6Q|) 
requires the construction of all possible parton shower histories [<&(£)] and the computa- 
tion of their probabilities "P(Ej). Since the number of parton shower histories £; for n final 
state particles grows like n\, this quickly becomes computationally expensive. 
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In general, we can assign 

d(£) = ^§-, (61) 

such that ^jQi(Si) = 1 automatically, and choose a(E) such that * s more readily 

calculable. The choices a(S) = 1 and a(S) = 7 ? (E)J(E) recover Eq. ( CT ) and Eq. (pD|), 
respectively. The better approximates J(S) the more uniform the total weights 

will be, resulting in a better reweighting efficiency At the same time, it is desirable to 
have a fast way to compute Yli a C^i)i such that the computational time used to correct the 
overcounting issue is less than (or at least comparable with) the computational time needed 
to evaluate the matrix element c($). In the next subsection we show that, with a suitable 



choice of a(£), the ALPHA algorithm of Refs. |77|, |78| can be applied to calculate the sum 
^Va^Ej), reducing the n\ problem to a 2 n problem, which is a significant gain. 

Before going on, we want to point out that the overcounting factor d(S) has appeared 
in other Monte Carlo programs, most notably MadEvent. Indeed, our solution to the phase 
space overcounting by the parton shower was inspired by the technique used in Ref. |L4| to 



slice phase space into complementary integration ranges. Given a set of n Feynman diagrams 



\M\- 



Em 



i=l 



(62) 



MadEvent performs independent phase space integrations over n different functions 

\M l \ 2 

fi = Wi\M\ 2 with Wi= — , (63) 



such that the sum over the n integrations ^Wil-M] 2 = \A4\ 2 gives the desired integrated 
squared- amplitude. The "overcounting factor" wi is not gauge invariant, but it does have 
the singularity structure of the z-th Feynman diagram, allowing a phase space integration 
that is aware of the relevant poles. GenEvA does essentially the same, but now a(E) will 
capture the singularity and color factor structure of a given parton shower history. 



B. The ALPHA Algorithm 

We want to calculate a quantity that corresponds to a sum over all possible graphs. This 
same issue arises in the computation of quantum amplitudes which are sums over all possible 
Feynman diagrams. Therefore, we can use known techniques to efficiently calculate the sum 
tt(Sj) by constructing a toy "Lagrangian" whose "generating functional" provides the 



desired "amplitude" as in Fig. [L5 



One simple algorithm for doing this is the ALPHA algorithm J77| , |75| , which is the core 



algorithm used in the ALPGEN event generator ||13|| . The ALPHA algorithm essentially works 
by reducing common subexpressions in the sum over all graphs. In a more familiar language, 
it exploits the fact that the calculation of a tree-level amplitude (which for us is the exact 
result) can be reduced to recursively solving equations of motion in the presence of external 
sources. 

Consider the following Lagrangian: 

£({<M) = 7^C ab (j) a (f) b + ^ D abc (j) a (j) b (j) c , (64) 
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FIG. 15: Schematic picture of the ALPHA algorithm. Given a quantity that can be described as 
the sum over all diagrams, one can create a Lagrangian C whose generating functional Z yields 
the appropriate sum over "Feynman diagrams". The ALPHA algorithm gives an efficient way of 
computing these sums as long as the fields (pi only appear linearly in C, hence the need for the off- 
diagonal propagator C~ b l = 5 a -b- In the above example, particle 1 could correspond to a photon, 
2 to a quark, 3 to a gluon, and 4 to an anti-quark. The 125 and 146 vertices could correspond to 
the electromagnetic couplings, the 34, —5 vertex to the q —* qg splitting function, and the 23, —6 
vertex to the q — > qg splitting function. 

where A is a dummy coupling constant and a are simple fields with no space-time depen- 
dence (i.e. field theory in + dimensions). The generating functional Z({J a }) of this 
theory is 6 

Z({ JJ) = -L({04) + <PaC ab J b , (65) 
where the fields a satisfy the variational condition (the equations of motion) 



5L 



The n-body correlation function is equal to 



CabJb 



>192 '-<Pn 



) = lim 



d n Z({J a }) 
Ji^o dJidJ 2 ■ ■ ■ dJ n 



(66) 



(67) 



and is identical to the sum over all Feynman diagrams (or in our case, the sum over all 
parton shower histories) with propagators C~ b l , vertices D a b c , and n external legs. 

To calculate Z at tree level, we can solve Eq. (|66| ) for the fields a . We start by expanding 
a in a series in A 



a = 0° + A0 a + --- + A 



n-2 in-2 



(68) 



6 The factor of C a b is included in the source term such that J a is a source for <fr a , not a source for C ob Vb- 
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where <p a is the coefficient of <p a at order \ % , and the series need only go to order X n ~ 2 , 
because an n-point correlator needs only n — 2 three-point interactions. The equation of 
motion following from Eq. (EBI) is 

<Pa = J a~ ^ C^D bcd (f) c (J) d , (69) 

which we can easily solve recursively order by order in A using Eq. (|68| ). The result is 

0a = ~\c- b l D bcd + 0^" 2 + ■■■ + 0*-V8) = ~\c- b l D bcd M , (70) 

i+j=k— 1 

which can then be plugged back into Eq. fl65l) . 

In the special case that £ is only linear in the <p a fields, the expression for the correlator 
in Eq. ( |67j ) simplifies. For n external fields (f>±, 02, • ■ ■ ,4> n - 

(01 • • • 4>n) = Z(J X = 1, J 2 = 1, • • • , J n = 1, Mother = 0) . (71) 

We will see below that the "Lagrangian" relevant for GenEvA is linear in the a , so Eq. ( |7T|) 
does indeed hold. 



C. A Shower "Lagrangian" 

To build a Lagrangian useful for the GenEvA algorithm, we introduce a separate field 
for every particle in a parton shower history (both for internal and external fields) and label 
each field by the specific four momentum p a and additional quantum numbers f a (flavor, 
spin, color, etc.) of the particle it represents, 

= 0/W a = 0a • (72) 

We stress that p a and f a are just labels; the field itself does not carry any momentum or 
quantum numbers. Note also that as far as the parton shower is concerned, all particles are 
distinguishable by their "L" and "i?" labels. In addition, we define all momenta as outgoing 
and introduce two separate fields for the same particle going into a vertex and coming out 
of a vertex, with the convention that if <p a represents an outgoing particle with outgoing 
momentum p a then 0_ a represents the same particle incoming with outgoing momentum 
— p a . The reason to have two fields for every particle is in order to use the result of Eq. (|7l|) 
which requires C to be linear in <fr a . 

Assume we have a parton shower history with n external (incoming or outgoing) particles 
with momenta Pi, where i runs from 1 to n. We again define all Pi as outgoing, such that 
overall momentum conservation reads Yli Pi = 0- For example, for a final state parton 
shower with n — 1 outgoing partons we include an "outgoing" center-of-mass particle with 
momentum 

n-l 

P n = -J2Pi. (73) 

i=l 
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The crucial point is that the momenta Pi span all possible momenta p a that can occur in 
any parton shower history with the same set of external particles. Explicitly, we introduce 
2 n — 2 different fields <p a corresponding to the momenta 

n 

p a = Y,j a Pi with 4 = 0orl, (74) 
i=i 

excluding the case where all c a are or 1, which would give p a = 0. The index a runs over 
1 < \a\ < 2 n ~ 1 — 1, and by definition we take c l _ a = 1 — c a such that 

n n 

P-a = J> - C\)P = - J2 4fi = ~P. , (75) 
i=l i=l 

as defined above. For simplicity, we also define c a = 5 % a for 1 < a < n such that the first n 
fields <f) a represent the external particles of the history. 7 Hence, while the number of parton 
shower histories grows like n\, the number of 4> a fields needed to calculate the sum over 
parton shower histories only grows like 2 n , and is therefore computationally tractable for 
reasonable numbers of external particles. 

To include the additional quantum numbers f a in the above discussion, we can replace 
Pa (p, f)a and Pi — > (P, f)i everywhere, supplementing Eq. ([F4]) with appropriate quantum 
number information (e.g. to include flavors). We only need to include a field <p a if there is at 
least one way to perform the specific sum over /. For example, if we have two outgoing quarks 
with momenta P\ and P2 then we do not need a field with momentum label pu = Pi + P2 
because there is no particle that can decay into two quarks. On the other hand, we can add 
two quarks and an antiquark of the same flavor by first adding the quark and the antiquark 
to give a gluon which can be added to the remaining quark. 8 

The addition of momenta is always unique. However, including additional quantum 
numbers, the final result of the sum over / can be ambiguous. For instance, including 
QED radiation and ignoring color, adding a quark and an antiquark can give a gluon or 
a photon. In this case, we either have to introduce a separate field <fi for each final result, 
or alternatively, include the additional quantum number, color in this case, that breaks the 
ambiguity. 

With the above conventions, momentum conservation in the parton shower requires 

Cab = C(a)5 a -b , C~t = cJ^jSa ,-b , Dabc = for p a + p b + p c ^ . (76) 

The first condition on the propagator means that an outgoing particle with p a can only 
connect to its incoming version with outgoing momentum p_ a = — p a , while the second 
condition enforces momentum conservation at each vertex in the parton shower history. 
Overall momentum conservation now implies that Z can only be linear in the sources for 
the external fields so Eq. ([FT]) holds. 9 Note that this relies crucially on the fact that we 



7 Except for the trivial case n = 2 where pi — P\ and p-i = P-2- 

8 In practice, the required fields 4> a are constructed recursively by always combining two previously con- 
structed fields, so adding flavors is trivial and requires no additional computational cost. 

9 To see this, assume that Z is quadratic or higher in one of the sources for the external fields with 
momentum Pi. This implies that there exists a correlator which has two or more copies of one of the 
original external fields. By construction, this cannot be the case, because it would require at least one 
field with momentum label containing 2Pi for some i. 
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introduced separate fields for incoming and outgoing particles, such that C ab would be off- 
diagonal. 

To obtain the desired sum over all parton shower histories we can now simply use Eq. ( |7T|) 
to plug the expansion of the fields in Eq. fl68|) back into Eq. fl65|) and extract the order X n ~ 2 
piece: 

^a(Ei) = (0i ■ ■■(/)„) 

1 1 ( 77 ) 

= "2 E C («)tt-g E ^010^ + C(a)0°0"; 2 , 

i+j=n— 2 i+j+k=n— 3 

where a sum over a, 6, c is understood and we also used Eq. ( [7(f ) for C a fe. The expressions 
for the fields 0^ in Eq. ( |77| ) are obtained from Eqs. (170) and Q75| ) and setting J^...^ = 1 and 
all other sources to zero, 



1 for 1 < a < n , 
otherwise , 



\—D^ bc Y, M- (78) 



2C(a, 

K ' i+j=k-l 

Once C(a) and D abc are specified, then Eq. ( [77] ) together with Eq. ( [78] ) is our final result 
for ^^(^j)- To obtain a(S) for a particular parton shower history (i.e. the numerator in 
Eq. fl6l|) ) we can simply multiply the appropriate factors of 1/C(a) and D abc for the given 
parton shower history. We will discuss the precise from of C a and the nonzero D abc next. 



D. Choice of Overcounting Function 

As discussed in Sec. |5 A] , ideally a(£) should equal the probability for generating a given 
history, but if we want to use the ALPHA algorithm outlined above, we have to be in a situation 
where that probability can be cast into a propagator C(a) _1 and a local vertex D abc , which 
can only depend on the momenta and quantum numbers associated with the indices a, b, c. 
As we have seen in Sec. ^ the splitting of a mother into two daughters requires information 
from the grandmother, which cannot be encoded in the local vertex D abc . Also, from the 
computational point of view, calculating the probability for a history requires numerically 
expensive Sudakov factors. Therefore, the compromise we will take is to choose 

C{a) = 1 (79) 

and have all nonzero D abc equal to a Jacobian-improved splitting function. 

We stress again that we are free to make any choice for C(a) and the nonzero D abc . 
For example, choosing C(a) = 1 and D abc = 1 the algorithm simply counts the number of 
parton shower histories, which is equivalent to using Eq. (|59|) and provides a cross-check on 
the implementation. The choice of overcounting "Lagrangian" only determines the efficiency 
of the program and not the resulting distributions. 

By definition, if D abc is nonzero then the flavors f a , f b , f c form a valid splitting with 
Pa + Pb + Pc = 0. In a parton shower, the mother particle is treated differently from the 
daughters, so to determine which of the three particles acts as the mother of the splitting, 
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we calculate the invariant masses t a ,b,c = V 2 a bc- Since the parton shower is ordered in t, we 
know that the largest of t a ,b,c belongs to the mother. Assuming for simplicity that this is t a , 
we use 

Dabc = ffa^hfc (ta, COS 9 a , (j) a ) J [t a ] t b , t c ) . (80) 

where J{t a ]tj,,t c ) is the Jacobian defined in Eq. (|)5|) and // a ->/ b / c is the splitting function 
defined in Section |3 B| . 10 Since this choice captures all the singularities of QCD, we can 
expect the corresponding to be reasonably close to "P(S)J(S). 

If t a = Eq M , then the splitting corresponds to the initial hard interaction, in which case 
we use 

D™ = J(t a ,t b ,t c ). (81) 

We can also add extra matrix element information to such as angular information or 
charges. In practice this is only important if there are two or more different hard-interaction 
vertices that can contribute to the same final state, such as e + e _ — > uudd, which can have 
both a e + e~ — > uu and e + e~ — > dd core. 



6. THE MATRIX ELEMENT 



The last missing ingredient in calculating the weight for a given parton shower history 
is the desired distribution cr [$(£)]. The appropriate choice for this distribution is entirely 
determined by physics considerations, and one should implement the best possible theoretical 
calculations available. Since this work focuses on the generation of the phase space, we will 
only consider two simple theoretical distributions as examples and discuss more elaborate 
choices in the companion paper 

We briefly discuss phase space issues for tree-level (LO) matrix elements cr( < I ) ). We then 
consider tree-level matrix elements improved to leading- logarithmic accuracy (LO/LL). This 
corresponds to the type of calculation one does in parton shower/matrix element merging |25|, 

|33| , [34], [35], |36|, [37|. To arrive at a simple expression for the LO/LL 
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cross section, we will use the notion of <Jj($) and the alternative form of the event weight 
in Eq. (ITT 



A. Tree-Level Matrix Elements 

The simplest theoretical distributions one could use in an event generator are given by 
fixed-order tree-level matrix elements. These results depend only on the scalar products 
of the four-momenta of the external particles and are defined entirely by Lorentz-invariant 
phase space variables $. While such tree- level calculations give at best a rudimentary 
description of the true QCD distributions, we will use them to validate our phase space 
generator in Sec. |7 A| and test its efficiency in Sec. 



10 There is a technical issue that the splitting functions are functions of It^ art , so strictly speaking there 
is not enough local information to even determine the correct splitting functions, because in the shower 
algorithm we set Eff rt to the maximum energy of the grandmother. However, regardless of which Eff Tt 
we choose to define D a b c , the singularities will still be correct. In practice, we calculate Eff 11 ' 1 assuming 
that the mother's sister is unbranched and that the grandmother's rest frame is the CM frame. 
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FIG. 16: 3-body phase space for e + e 



qqg. Region 2 corresponds to the 2- jet region because 

However, the 



the parton shower will not populate regions with both tqg and tqg less than tcut 
parton shower can "accidentally" distribute events in the 2 regions where either t qg or tq g is less 
than t cut . These regions must be vetoed because they occupy a QCD singularity away from a 
shower singularity. While there is no QCD singularity associated with t q q < t cut , it is convenient to 
also veto that region for simplicity (e.g. when comparing to MadEvent in Sec. [7 A| ). The remaining 
part of phase space defines the 3-jet region 3. 



We want to touch on two additional subtleties having to do with phase space vetoes. 
Recall that at tree level, there are singularities in QCD matrix elements when two partons 
get close in Uj = (pi+pj) 2 . In full QCD, these singularities are regulated by Sudakov factors, 
but a tree-level event generator has to cut off these ty singularities by hand. In the context 
of the GenEvA algorithm, the natural phase space cut variable is virtuality. After all, if we 
run the shower down to t cut , then the shower never encounters the singular region of QCD. 

However, when we run the shower, it is possible for two particles to "accidentally" have 
virtualities smaller than t cut , and the matrix element will return a singular answer. This 
sickness also appears in the overcounting factor, because the ALPHA algorithm will encounter 
a value of D abc that cannot be generated by a shower with a t cut ending scale. Consider 
3-body phase space for e + e~ — > qqg as in Fig. 16. There are regions of phase space for which 
a large angle emission from q — > qg can put the gluon closer to the anti-quark than t cnt . The 
q — > qg splitting cannot generate those same four-momenta because it is below the cutoff 
scale, so there is no "partner" shower history for the event. Moreover, the shower generated 
those four-momenta far away from any QCD singularity, while the four-momenta are indeed 
in the singular region, so the event weight will be very large. 

The solution, of course, is to veto events from the shower for which the minimum tij 
between singularity-producing partons is less than t cut . Note that we need not veto events 
for which the q and q are closer than t cut , because in the context of e + e~ — *> qqg, there is 
no qq singularity. 11 If for simplicity we do not wish to impose a flavor-dependent t cn t scale, 
then we can always veto events where any pair of particles (singular or not) are closer than 



11 There is a singularity in e + e — * qqq'q' , of course. 
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t cut . Denoting the fixed-order tree-level matrix elements by cr tree ($), this is equivalent to 
taking 

a LO ($)=a trcc ($) l[9(t ij >t cut ). (82) 

u 

Using Eq. (|12"D, the corresponding event weight is 



which is what we will use in Sec. 17 A. 



The second subtlety has to do with the fact that the scale t matc h is ambiguous for n max - 
body phase space. Recall from Sec. |3 Ej that to create fully exclusive events with no dead 
zones, we want to start a parton shower approximation at t ma tch- When n max -body phase 
space is generated, the scale t ma tch will vary on an event-by-event basis to ensure that the 
attached parton shower will populate all of {n > n max )-body phase space. However, since 
the same n max four-momenta can be created by different shower histories during the phase 
space generation, they can be associated with several different t ma t c h scales. 

As far as phase space vetoes are concerned, the variable t mat ch can simply be ignored, 
because even after truncating to {t cut , n max }, the shower does populate all of n max -body 
phase space down to t cut . However, if we start a parton shower approximation at t ma t c h, we 
do have to be mindful that events having the same n max partonic four-momenta will start 
the shower at different scales. Far from being a bug, this t mat ch ambiguity will be exploited 
in the next subsection to create the leading-logarithmically improved event sample. As we 
discuss in the companion paper [[|, this ambiguity is a physical ambiguity that has to be 
resolved with calculations. 



B. Logarithmically-Improved Matrix Elements 

It is well known that fixed-order calculations do not give an adequate description of true 
QCD distributions in singular regions of phase space. This is because fixed-order calculations 
only reproduce the power singularities that occur when intermediate propagators go on shell. 
However, in addition to these power singularities there are logarithmically enhanced terms, 
which start to show up one higher order in perturbation theory. These logarithmically 
enhanced terms have the form [a s log 2 (t/£ , Q M )] n . Thus, the effective expansion parameter 
in the singular regions of phase space is not a s , but rather a s log 2 (t / E"^ M ) . These double- 
logarithmic terms can be resummed to all orders in perturbation theory, giving an effective 
expansion parameter a s \og(t/Ec M ), which is only enhanced by a single logarithm. However, 
such resummed calculations are more difficult to implement, since they not only depend on 
the four-momenta of the external particles, but also on a choice of scales which determine 
the exact form of the resummation. 

Resummed calculations are crucial if one wants to merge partonic calculations with parton 
shower algorithms to produce fully exclusive events, as required for meaningful comparisons 
with experimental collider data. The parton shower automatically includes the resummation 
of the double- logarithmic terms via the Sudakov factor discussed in Eq. (P§D, and therefore 
has double-logarithmic sensitivity to the scale at which it is started. Thus, the absence of 
resummed partonic calculations leads to a double-logarithmic dependence on the unphysical 
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matching scale t m atch between the partonic calculation and the parton shower. This issue is 
discussed in much more detail in the companion paper [|TJ. 

Most implementations of resummed calculations in event generators define the relevant 
scales in terms of the external four-momenta. An example is the so-called CKKW prescrip- 
tion to merge fixed order calculations with parton showers. In this prescription, one 
first generates an event based on tree-level matrix elements and then uses a flavor-aware hp 
algorithm on the external particles to find the most likely parton shower history that could 
have generated this event. The resummed partonic result is then obtained by reweighting the 
fixed-order event by the appropriate no-branching probabilities associated with this specific 
parton shower history. 12 

There are two main issues with this approach and other similar ones. First, this kind 
of algorithm only allows partonic calculations which are the simple product of a tree-level 
matrix element and Sudakov factors. In general, the resummed expression can contain vari- 



ous terms with differing logarithmic structures. A recent study |7(J, |7T| of the production of 
highly energetic partons in soft-collinear effective theory (SCET) p6[ |68|, |69| has shown 
that the logarithmic structure is obtained by a sequence of matching and running calcula- 
tions, which leads to a sum over terms, each with a different resummation of logarithms. 
These different terms can be identified with the different possible shower histories, which is 
due to the fact that the effective theory and the parton shower are both derived from the 
soft-collinear limit of QCD. 

Second, reweighting tree-level calculations by Sudakov factors may result in inefficient 
event generators. This is because the product of many Sudakov factors can give numerically 
very different results for different events depending on the precise scales found by the kp 
algorithm or any other jet identification procedure. Thus, if the initial set of events was 
distributed only with the correct power-singularity structure, the resulting events after in- 
cluding logarithmic information can have widely differing event weights, leading to a poor 
statistical efficiency. Of course, one can easily avoid this issue by accounting for both the 
power and logarithmic singularities at the time of initial event generation. 

In GenEvA both of these issues are solved automatically, because the underlying phase 
space generator is itself a parton shower. First, we can assign different weights for different 
parton shower histories, allowing an easy implementation of the sum over different logarith- 
mically resummed terms. Second, the events are directly generated based on the singularity 
structure of QCD and already include the leading-logarithmic effects. Thus, reweighting 
the parton shower to resummed partonic expressions will actually give much more uniform 
weights than reweighting to fixed-order matrix elements. We will come back to this issue in 
Sec. 7~C when discussing the efficiency of GenEvA and focus for now on how to implement 



resummed partonic calculations in GenEvA. 

A particularly useful expression for resummed partonic distributions is given in the com- 
panion paper ffl]. Starting with Eq. (0), note that the parton shower by itself is creating 
events with w(H) = 1. Therefore, the internal GenEvA parton shower distributes events with 
an effective ^-dependent "matrix element" 

^shower = J( S ) = Q(S)A(E) , (84) 

where Q(H) is the product of a hard scattering matrix element, splitting functions, and 



12 CKKW also corrects for the running of a s . 
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Jacobians, and A(E) is the product of Sudakov factors. Starting from the fixed-order tree- 
level matrix elements cr tree ($), an LO/LL sample can be defined by (see Eq. (Ill) in Ref. [|TJ) 

a, LO/LL [$(S)] = cr tree [$(E)] - A(S) , (85) 

where the sum runs over all shower histories £«($) that yield the same phase space point 
$. Therefore, the formula for the LO/LL event weight is quite simple 13 

w lo /LL(s) _ ^ LO/LL [^)] _ ^*<m m 

W V(E)J(E) -XiQPiWV ( } 

Although Eq. (|85D includes a Sudakov factor, the LO/LL event weight no longer has any 
Sudakov information in it, because the internal GenEvA shower already distributes events 
with the correct double-logarithmic information. Since the Sudakov factors A(£) contained 
in P(E) in the LO weight in Eq. fl83|) come at some computational cost, creating this LO/LL 
merged sample is computationally even more efficient than the LO result alone. 

Note the difference between this result and CKKW. Here, we are using every shower 
history to define an appropriate Sudakov factor. The relative contribution of each Sudakov 
factor is given by 

gg) (* 7] 

E,Q[Ei(*)]" 1 } 

In those regions of phase space that are dominated by one shower history, one Q(Y) — > 
oo and therefore one Sudakov dominates. This corresponds to the approach of Ref. fl26| . 
In CKKW [f25|j , the chosen Sudakov factor always corresponds to the dominant history, 
whereas in Eq. (|85|) , away from the singular regions, an average Sudakov is used, more 
closely mimicking the SCET calculation fFO], [7T]| . 

It would be interesting to study which approach better corresponds to experimental data, 
though we emphasize that the GenEvA algorithm can of course distribute events that are 
morally equivalent to the CKKW procedure, by simply choosing the partonic cross section 

a CKKW ($) = a tree ($) A [£ doln (<t>)] , (88) 

where Sd m( < ^ ) ) is the dominant shower history for the phase space point <3>. 



C. Further Logarithmic Possibilities 

GenEvA also allows a numerical treatment of an interesting single-logarithmic issue. In 
the LO matrix element Eq. fl82|), we vetoed events for which two partons were "accidentally" 
closer than t cut . This made sense, because we were not using the Sudakov information in 



Since we use the ALPHA algorithm as described in Sec. |5| to compute the sum over parton shower histories, 
the precise splitting functions in Eq. ( |S5| ) are those used in Eq. (|8(]). The actual splitting functions 

used in the shower slightly differ from these by the choice of E max . Therefore, in practice, Eq. ( |86| ) 
contains an additional factor 
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the LO calculation. However, imagine the case that the tree-level matrix element were 
identically equal to the product of splitting functions 



( g g) 

i 

such that the event weight in Eq. fl86|) is always one. Because the parton shower preserves 
probability, the < t cut veto will now decrease the total cross section for e + e~ — > n jets 
from the Born value. Putting in the correct expression for cr trco ($), this will result in a 
single-logarithmic decrease in the cross section. While this effect is formally beyond LO/LL 
order, it is numerically important for high-multiplicity matrix elements. 

As we argued in Fig. [16|, events that violate the i cut constraint should be thought of as 
actually contributing to lower-dimensional phase space. Imagine we take an event sample 
and truncate to {t cu t, ^max}- As explained earlier, this gives (n < n max )-parton events with 
^match = tcut, and n max -parton events for which t ma t c h changes from event to event. Even 
after truncation, there will be some n-parton events with 3 < n < n max left that violate 
the t cnt constraint. Instead of vetoing these "rogue" events, we can further truncate them 
to some scale t' > t ma tch until the t cut violation disappears. This will take all rogue n- 
parton events and make them into (< n)-parton events, recovering the cross section that 
would otherwise be lost to the t cut veto. We will refer this additional truncation up to t' as 
"truncation-prime" . 

Remember that these events were originally truncated to {t cut , n max }, so we have to restart 
the shower at the originally determined £ m atch scale (and not t') , otherwise we would double- 
count existing regions of phase space from the ordinary n-parton events. If the probability of 
getting a given n-body shower history E n after truncation is P trunc (S n ), then the probability 

ptrunc ( Sn ) tQ get the 

same shower history after truncation-prime is given schematically by 

p tmnc '(s n ) = p trunc (s„) + d ( s «W£n) P trunc (s n+J ) , (90) 

i=l 

where d(S n+ j/S n ) represents an integral over those (n + i)-body shower variables that are 
not shared by the n-body shower, and n represents those regions of (n-M)-body phase space 



that really look like n-jet events (e.g. the 2 regions in Fig. [If]). Note that for n = n max , 
truncation-prime does not affect the event probability. 

If we now use the same event weight as in Eq. (^6|) for the resulting event sample after 
truncation-prime, then this is equivalent to choosing 



^ [4(E)H ^ [4(E)] _|^_ A(E) ^^. (91) 



In other words, we can do the integrations in Eq. ( P0{ ) numerically by truncating the shower 
to t'. This p trunc /ptrunc en h ancemen t is formally beyond leading- logarithmic order, so it 
is safe to include it. Numerically, it helps yield more stable total cross sections as t cut is 
varied. The main downside to using the numerical 7 ?trunc /-ptrunc me thod is that it breaks the 
reversibility of further truncation. That is, once truncation-prime is applied, the probability 
for obtaining a set of four-momenta is proportional to p trunc and thus extremely convoluted, 



so one can no longer use the argument at the end of Sec. |3 E| that reweighting can occur 
after hadronization. 
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7. RESULTS 



In this section, we show results for the process e + e~ — > n jets using the partonic cross sec- 
tions defined in the previous section. As input for the fixed-order tree-level result <7 tree ($), we 
use MadGraph [79] to obtain numerical Fortran HELAS routines for e + e~ — > n partons 14 , 



where 2 < n < 6, and a "parton" is either a gluon or a u, d, s, or c quark (neglecting quark 
masses). In this section, we always use the center-of-mass energy and phase space cut 

E CM = 1000 GeV , = 100 GeV , (92) 

unless otherwise noted. 

We first validate the GenEvA algorithm by using the simple LO cross section from Eq. ( |52"D 
and comparing the output of GenEvA with that of MadEvent [TH] . We then show some rudi- 
mentary results using the LO/LL cross sections in Eqs. (^) and (|9"TD, with a more thorough 
discussion given in the companion paper [p]]. Finally, MadEvent offers a useful benchmark to 
compare the efficiency of the GenEvA algorithm. We will see that GenEvA's efficiency is at or 
above MadEvent levels, and that GenEvA itself is more efficient at distributing logarithmically- 
improved LO/LL results than pure tree-level LO results. 



A. Comparison with MadEvent 

Since GenEvA uses the same MadGraph/HELAS matrix-element engine as MadEvent, both 
programs should and do give identical results. For the purpose of this comparison we are not 
attaching any additional parton shower emissions, and simply compare the partonic cross 
sections for a fixed number of final state partons. 

We first consider the results for the total cross sections, which are shown in Table [I 
for the n-parton samples combining all flavor subprocesses as well as for various individual 
flavor subprocesses. One can clearly see that the two programs agree within statistical 
uncertainties. Notice that while in MadEvent the relative uncertainty for each subprocess 
is roughly the same, in GenEvA the relative uncertainty scales as the square root of the 
cross section for the given subprocess. This is because GenEvA automatically produces a 
sample containing the sum of all subprocesses in a single run, with the ratio of the various 
subprocesses roughly given by the QCD symmetry structures. This means that the most 
important channels automatically get the most statistics and thus the smallest relative 
uncertainties, and no extra integration time is spent on less important channels. In MadEvent 
the individual subprocesses are generated separately and combined afterwards. 

Next, we compare the differential distributions generated by both programs for a few 
simple partonic observables. We have checked a large range of partonic observables and 
in all cases the distributions obtained from GenEvA and MadEvent agree perfectly within 
the statistical uncertainties, verifying that the GenEvA phase space generator is working 
properly. As an illustration, we show in Fig. O the distribution of the pairwise invariant 



mass between the 2nd and 3rd hardest partons for 3- and 4-parton matrix elements (left 
panel) and between the 4th and 5th hardest partons for 5- and 6-parton matrix elements 



14 For simplicity, we are only including the diagrams with e + e — > 7* — > partons. An intermediate Z-boson 
could easily be included. 
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process 


MadEvent 


GenEvA 


process 


MadEvent 


GenEvA 


LO 3 (fb) 


216.71 ±0.21 


216.77 ±0.22 


LO 5 (ab) 


2542 ± 3 


2543 ± 3 


uug 
ddg 
ssg 
ccg 


86.62 ±0.13 
21.75 ±0.07 

21.63 ±0.06 
86.71 ±0.13 


86.60 ±0.18 
21.55 ±0.10 
21.73 ±0.10 
86.70 ±0.18 


uuggg 

ddggg 
uuddg 
uuuug 


912 ±2 
227.5 ±0.9 
33.8 ± 0.2 
25.6 ±0.2 


912 ±2 
228.3 ± 0.8 
34.3 ± 0.4 
25.7 ±0.3 


LO 4 (fb) 


36.44 ±0.04 


36.49 ±0.04 


LO 6 (ab) 


67.9 ± 0.3 


68.0 ± 0.2 


uugg 
ddgg 
uudd 
uuuu 


14.00 ±0.03 
3.504 ±0.013 
0.175 ±0.001 
0.132 ±0.001 


14.00 ±0.02 
3.511 ±0.011 
0.180 ±0.003 
0.132 ±0.002 


uugggg 
uuuugg 
uuuuuu 
uuddss 


22.41 ±0.09 
1.117 ±0.006 
0.005 ±0.001- 
0.019 ±0.001" 


22.29 ±0.12 
1.14 ±0.03 
0.005 ± 0.001 
0.020 ± 0.005 



TABLE I: Comparison of total cross sections between GenEvA and MadEvent. Both programs agree 
within statistical uncertainties (0.001 - indicates an error too small to be reported by MadEvent). 
In MadEvent, the relative uncertainty on each sample is roughly constant, whereas in GenEvA, less 
populated channels are allowed to have higher relative uncertainty. This happens because the 
ratio of various subprocesses is given in GenEvA by the QCD symmetry structures, and no extra 
integration time is spent determining the ratio between different channels. 



(right panel). Here, hardest is denned by the parton with the largest energy. We will discuss 



the different efficiencies of the two programs later on in Sec. 7C. 



B. Logarithmically-Improved Results 

The main usefulness of an LO/LL partonic calculation is for creating fully merged LO/LL 
samples interfaced with a parton shower. We do this in the companion paper fl|] and here 
only give partonic results, which by themselves have only limited usefulness, but show the 
effect of Sudakov suppression. 

From Eq. (|85|) , we see that the LO/LL partonic calculation depends on a Sudakov factor, 
which in turn depends on the t ma tch scale of the event after truncation. For an n-body matrix 
element with fixed n, there are two different cases, as discussed in Sec. |T|: If n = n max , 
then tmatch > t cut is an event-specific scale corresponding to the last emission and the result 
should be regarded as an inclusive sample. If n < n max , we always have t ma tch = t CVLt and the 
result should be regarded as exclusive. For n < n max , we also have the choice of whether 
to include the truncation-prime correction factor in Eq. (|91|). These three possibilities are 
labeled as 15 

LO/LL inc. n = Eq. (|85|) with n = n max , 
LO/LL exc. n = Eq. (|85|) with n < n max , 

LO/LL' exc. n = Eq. (|l]) with n < n max . (93) 



15 In the companion paper Q, we only show LO/LL inc. and LO/LL' exc. results and will therefore remove 
the prime symbol from the LO/LL notation. 
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FIG. 17: Comparison between MadEvent and GenEvA for 3-, 4-, 5-, and 6-parton tree-level matrix 
elements. Left panel: The pairwise invariant mass of the 2nd and 3rd hardest partons for 3- 
and 4-parton matrix elements. Right panel: The pairwise invariant mass of the 4th and 5th 
hardest partons for 5- and 6-parton matrix elements. The two programs agree within the statistical 
uncertainties, validating the GenEvA algorithm since both programs use the same HELAS/MadGraph 
matrix elements. Note the hard cutoff at y/t cnt = 100 GeV. 



Differential distributions for these three possibilities are shown in Fig. for n = 3, i.e. 
e + e~ — > qqg. Again, the hardest parton is defined as the parton with the largest energy. 
When the invariant mass between the 1st and 2nd hardest partons is at the lower threshold, 
the three partons are well separated, with a large emission scale for the gluon. Since the 
LO/LL inc. sample only captures the logarithms between Ecm and the gluon emission scale, 
it contains no Sudakov suppression in this part of phase space, and hence asymptotes to 
the LO result. The LO/LL exc. sample has an additional overall Sudakov suppression from 
leading-logarithmic running from the gluon emission scale down to the fixed matching scale 
^match = t cnt - The upper endpoint is the 2-jet region with a low gluon-emission scale, so all 
three LO/LL samples are Sudakov suppressed compared to the LO sample. The LO/LL' 
exc. sample has roughly the same shape as the LO/LL exc. sample but the overall scale 
is increased, because the n integrations it contains restore some of the lost cross section 
from phase space vetoes. The reason it asymptotes to the LO/LL inc. sample at the upper 
endpoint is because in this region it is most likely for two partons to be "accidentally" 
closer than t cnt - The identical effect can be seen in the right panel when the invariant mass 
between the 1st and 3rd hardest partons is small. Note that the Sudakov suppression from 
the LO result to the LO/LL results is formally double- logarithmic, while the increase from 
the LO/LL exc. to the LO/LL' exc. result due to the inclusion of the n integrations is 
formally single-logarithmic. 



C. The Efficiency of GenEvA 



We argued in Sec. |1 Q that GenEvA has several advantages over more traditional algo- 
rithms, so it is interesting to see how its efficiency compares to other event generators. 
We begin by addressing the speed of GenEvA for fixed-order calculations by comparing it 
to MadEvent JTJj], and then show how GenEvA is optimized to generate resummed partonic 
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FIG. 18: The effect of Sudakov factors on 3-parton matrix elements e + e~ — * qqg. In the LO/LL 
inc. sample, the Sudakov suppression only captures the logarithms between Ecm and the gluon 
emission scale. In the LO/LL exc. sample, there are additional no-branching probabilities from the 
gluon emission scale down to the fixed t m atch = icut scale, giving effectively an "exclusive" 3-jet 
result with an overall Sudakov suppression. The LO/LL' exc. sample includes the n correction 
factor from Eq. (^), while the LO/LL exc. sample does not. Left panel: Invariant mass between 
1st and 2nd hardest partons. The LO and LO/LL inc. samples agree near the lower endpoint, 
because this is the 3-jet region where Sudakov suppression is not active. The upper endpoint is the 
2-jet region, and the n integrations in the LO/LL' exc. sample restore lost cross section compared 
to the LO/LL exc. sample, asymptoting to the LO/LL inc. result. Right panel: Invariant mass 
between 1st and 3rd hardest partons, which shows similar effects. In the companion paper [Q], we 
will combine the LO/LL partonic results to create a merged LO/LL sample. 



results. 

Comparing these two programs is not straightforward, since the different methods used for 
phase space generation make it difficult to design a fair comparison. First, GenEvA is written 
in C++, whereas MadEvent is written in Fortran. Thus different levels of optimization in the 
compilers can affect the results. However, given that both programs use the same Fortran 
MadGraph/HELAS matrix elements and because evaluating the matrix element is the slowest 
part of GenEvA, the comparison is reasonable though not optimal. Second, in a single run 
GenEvA directly produces a fully inclusive sample, containing a sum over all subprocesses 
and multiplicities, whereas MadEvent generates events for each subprocess and multiplicity 
separately and combines them at the end of the run. Thus, MadEvent is producing many 
events that are thrown out at the end of the day. This is a source of inefficiency that could 
easily be ameliorated by first calculating the approximate cross sections for the various 
subprocesses and then generating only the events required. In our speed comparisons, we 
therefore give timing information for n-parton samples with fixed n and summed over all 
subprocesses for both GenEvA 16 and MadEvent, as well as for MadEvent samples with only 
one flavor subprocess. 



To obtain an n-parton sample with fixed n in GenEvA we simply set n — n max and discard events with 
fewer than n final state partons. 
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T eff (msec) 


T .9 (msec) 


GenEvA LO 3 




0.789 


0.57 


0.62 


GenEvA LO/LL inc. 


3 


0.965 


0.47 


< 0.47 


MadEvent 3 




0.982 


2.6 


< 2.6 


MadEvent uug 




0.994 


3.0 


< 3.0 


GenEvA LO 4 




0.525 


1.7 


2.2 


GenEvA LO/LL inc. 


4 


0.713 


1.3 


1.5 


MadEvent 4 




0.809 


11.1 


11.4 


MadEvent uugg 




0.752 


5.4 


5.7 


GenEvA LO 5 




0.390 


10.0 


15 


GenEvA LO/LL inc. 


5 


0.557 


8.6 


10.8 


MadEvent 5 




0.843 


62 


64 


MadEvent uuggg 




0.833 


27 


27 


GenEvA LO 6 




0.298 


160 


250 


GenEvA LO/LL inc. 


6 


0.396 


150 


230 


MadEvent 6 




0.809 


1900 


2300 


MadEvent uugggg 




0.784 


330 


350 



TABLE II: Comparison of the speed and efficiency between GenEvA and MadEvent. rj c ^ is the 
statistical efficiency defined in Eq. (|l7|), and T e s is the time required to create one statistical event. 
To.g is the time to create one statistical event in a partially unweighted sample with rj e s(wo) = 0.9, 
as defined in Eq. (|l9|). Note that GenEvA performs better on Sudakov-improved (LO/LL) results 
than on tree-level (LO) matrix elements. GenEvA is in all cases faster than MadEvent, in most 
channels by a factor of a few, despite the fact that GenEvA's efficiencies are lower since it is a 
fixed-grid algorithm. 

The efficiency results are shown in Table |J, where 3-, 4-, 5-, and 6-parton matrix elements 
are compared. Shown are the statistical efficiencies r] eS defined in Eq. (|i~7D, as well as the 
time T e fj needed to create one statistically equivalent event. Because GenEvA is a fixed- 
grid algorithm, its statistical efficiency r] e ff is lower than that of MadEvent, where grid 
optimization is used. However, the time needed to create a statistically equivalent weighted 
event sample is smaller for GenEvA than MadEvent, presumably because of the speed of 
the parton-shower-based phase space algorithm. To assess this efficiency/speed tradeoff, we 
also report T g, the time needed to create one statistically equivalent event in a partially 
unweighted event sample with r] eS (w ) = 0.9 as defined in Eq. (|19|). That is, we unweight 
just those events with weights below a threshhold weight wq, with wq chosen such that the 
final event sample has a 90% statistical efficiency. By this To.g measure, GenEvA is still faster 
than MadEvent in all channels, often by a factor of a few. 

We do not show comparisons of the unweighting efficiencies between GenEvA and 
MadEvent, because there is no reasonable way to do so. As argued in Sec. |2 C| , unweighting 
efficiencies are sensitive to the exponential tails of the weight distribution. Depending on 
how long one runs the two programs, one may or may not sample these tails, giving wildly 
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FIG. 19: Comparison of the weight distributions between GenEvA LO and LO/LL samples for 3- 
to 6-parton matrix elements. Because the GenEvA phase space generator already includes leading- 
logarithmic information, the LO/LL weight distributions are more strongly peaked at their average 
weights. More importantly, there is a considerable suppression of the high-weight tail when going 
from the LO to LO/LL samples. As expected, the effect of leading logarithms is less relevant 
for the higher- multiplicity samples, as phase space suppression keeps the double- logarithms from 
growing too large. 



varying unweighting efficiencies during one run and between different runs. Furthermore, it 
is completely sufficient for all practical purposes to have partially unweighted samples with 
a reasonably high statistical efficiency. 

As we discussed above, the fact that GenEvA is using a parton shower as a phase space gen- 
erator makes it ideal to distribute events according to logarithmically-improved partonic cal- 
culations. This is because the parton shower already contains the correct double-logarithmic 
behavior, such that the resulting weights are expected to be uniform, giving rise to a high 
statistical efficiency. To see this, consider the event weight given in Eq. fl8S|). Since the split- 
ting functions reproduce the tree-level matrix elements of QCD in the collinear limit, the 
denominator has exactly the same singular behavior as the numerator, with no logarithmic 
differences. 

This increase in efficiency for LO/LL event generation can be seen in a few different ways, 
the simplest being that in Table [TTJ the three different efficiency/ speed measures are all better 
for the GenEvA LO/LL samples than for the GenEvA LO samples. In Fig. [19|, we show the 
distribution of event weights between the LO/LL and LO samples. For n = 3 especially, 
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FIG. 20: Comparison of the statistical efficiencies rj e g between GenEvA LO and LO/LL samples for 
3- to 6-parton matrix elements. As the matching scale /x = \J t ma tch is lowered, the efficiency of the 
LO/LL sample stays roughly constant, while the efficiency of the LO sample drops dramatically. 



the weight distribution is far more peaked near the average weight in the LO/LL sample 
compared to the LO sample, showing that the inclusion of leading-logarithmic information 
softens the high- weight tails. Even for n = 6, where Sudakov logarithms play a much smaller 
role in determining the differential distributions, there is still a factor of 5 decrease in the 
number of high-weight events. The change in efficiency is shown even more dramatically in 
Fig. which compares the statistical efficiency between the LO/LL and LO samples as 
the matching scale \x is changed. For the LO/LL samples, the efficiency is nearly constant 
as jj, varies, while the LO samples exhibit a sharp decrease in efficiency for fi — ► 0. 

Thus, GenEvA is faster at distributing resummed results than tree-level results. We note 
that the efficiency of GenEvA could improve further with the use of a more realistic internal 
parton shower. For example, tree-level matrix elements include the effect of color coherence 
which is absent in the virtuality-ordered shower currently used in GenEvA. If a p^-ordered 
shower could be implemented as discussed in Sec. |S Q , then the leading color coherence 
effects would already be captured by the phase space generator, which should lead to more 
uniform weights. 
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8. CONCLUSIONS 



We have shown how to construct a phase space generator from a parton shower. The 
resulting GenEvA algorithm uses a very different strategy compared to more conventional 
approaches. Most available generators use a non-recursive, multi-channel, adaptive-grid 
algorithm to distribute events in phase space directly according to a desired distribution, 
which are often tree-level matrix elements. In contrast, GenEvA uses a parton shower to 
generate events, and in a second step reweights the events to the desired distribution. 

Though the GenEvA algorithm uses a fixed grid with no adaptation, the recursive nature 
of the parton shower allows for an incredibly fast generation of phase space points. Since the 
parton shower has the same singularity structure as QCD, the effectively fixed grid used in 
the GenEvA algorithm is presumably very close to the final adapted grid in traditional phase 
space generators. As a result, the GenEvA algorithm is able to compete with and likely even 
outperform traditional generators. 

The main advantage of the GenEvA algorithm arises when one goes beyond fixed-order 
calculations. It is well known that leading-logarithmic terms have to be taken into account 
in order to obtain realistic QCD predictions. In particular, when merging matrix element 
calculations with subsequent parton shower evolution, leading-logarithmic information is 
essential to avoid strong dependence on the unphysical matching scale. Unfortunately, ana- 
lytic expressions for logarithmically improved, fully differential distributions are prohibitively 
difficult to obtain for a large number of external particles, since the required scales in the 
logarithms depend in a complicated way on all kinematic invariants in the event. 

For this reason, traditional generators often implement leading-logarithmic improvements 
after the initial phase space generation is complete, and use the known kinematics of the 
event in order to calculate the required scales in the logarithms. In contrast, the parton 
shower used to generate phase space points in GenEvA automatically contains the leading 
logarithmic behavior of QCD by construction, and for this reason GenEvA is in fact more 
efficient in distributing logarithmically improved distributions than fixed-order calculations 
and thus has the potential to outperform traditional techniques. 

To conclude, we discuss some potential extensions to the GenEvA algorithm to make it 
useful for LHC physics. 

A. Towards Hadronic Collisions 

Currently, the GenEvA algorithm is only designed to work for e + e~ collisions with no 
initial state photon radiation. The technical reason for building GenEvA for e + e~ first is 
that the final state parton shower with fixed center-of-mass energy is easier to understand 
and verify. However, the name GenEvA is obviously inspired by the site of the upcoming LHC 
experiment at CERN. We anticipate extending the GenEvA algorithm to hadronic collisions, 
and one might wonder whether there are any additional technical hurdles to overcome. 

To build GenEvA for e + e~ collisions, there were three technical problems that needed 
to be solved. First, in order to generate events that conserve energy-momentum, parton 
showers need to adjust the kinematics of the parton splittings after they are generated, which 
typically makes it impossible to calculate the probability analytically. Second, different 
parton shower histories can give rise to the same points in Lorentz-invariant phase space, 
such that a parton shower covers phase space multiple times. Third, one needs a way to 
"truncate" a high-multiplicity event generated by the parton shower to a lower-multiplicity 
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event for which the probability is still known analytically. 

In this paper we explained in detail how these technical issues are solved in the GenEvA 
algorithm. The first issue is solved by employing the analytic parton shower developed in 



Ref. [|5!J. The second issue is solved by introducing the overcounting factor which multi- 
plies the generated event weight. Finally, the double-branch probability used in this paper 
enforces a global evolution, which allows for a simple truncation algorithm. 

These technical solutions should have straightforward generalizations to hadronic colli- 
sions. While we know of no analytically calculable initial state shower, building one should 
not be too difficult as it is well known that phase space can be recursively described both 
in terms of an s-channel "final state" tree and a t-channel "initial state" tree. The trick 
to creating the analytic final state shower was to use cos 9 as the primary variable instead 
of z, and an analogous trick should be possible for the initial state shower. The issue of 
double-counting will require figuring out how to deal with interference between initial state 
and final state shower configurations, but as a last resort, one can always do a brute force 
sum over parton shower histories to solve the double- counting issue. Truncation should be 
more-or-less identical to the final state shower, because many initial state showers are based 
on backwards evolution, and therefore can have a notion of well-ordering. 

Hence, the main new issue in hadronic collisions is the proper treatment of parton dis- 
tribution functions (PDFs). However, apart from issues of efficiency, this is really not an 
issue, because parton distribution functions simply introduce two new variables x\ and x 2 
and define a new choice of partonic calculation 

da = /xOci, Q?)/ a (z 2 , Ql) dx! dx 2 da 12 ^ x = </($') d$' , (94) 

where $' includes x% t 2 in addition to the ordinary phase space variables, and <j'(&) is the 
fully differential cross section including PDF information. The point of reweighting is that 
no matter how events are selected in events can be distributed according to &'(<&'), and 
it is the user's responsibility to, for example, use NLO PDFs with NLO calculations, use 
the proper scales for Q\ 2 , and so on. 



B. Heavy Resonances 

At the LHC, the most interesting Standard Model backgrounds to Beyond the Standard 
Model signals come not from QCD itself but from heavy resonances like W, Z, and t created 
in association with QCD radiation. Because GenEvA is based on a parton shower, one might 
wonder whether it is possible to include heavy resonances in the GenEvA algorithm. 

For this purpose, it is useful to think of GenEvA not as a parton shower but as simply 
a recursive phase space generator, and resonances corresponds to just a funny choice of 
"splitting function" . Because events will be reweighted at the end of the day to the desired 
distribution, all what is important is that the "splitting functions" have the right singularity 
structure, which for the purposes of heavy resonances is the Breit-Wigner line shape. 

What will be particularly bizarre about the GenEvA algorithm with heavy resonances is 
that there will be "Sudakov factors" associated with resonances, but as long as the resonances 
are narrow enough, this Sudakov factor will either be unity if t cnt > m 2 or zero if t cnt < tit 2 , 
meaning that a resonance is either on shell or it is decayed. This also means that even 
if t cu t > m? and a resonance is treated as on shell for the matrix element calculation, if 
the GenEvA shower is restarted at t cnt , the Breit-Wigner distribution will be recovered for 
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the resonance with the rest of the shower history compensating for the change in energy- 
momentum. 

Once the logistics for handling mass thresholds with heavy resonance is established, the 
treatment of not-so-heavy resonances like the bottom quark will be straightforward, and one 
would use the ordinary massive Altarelli-Parisi splitting functions for b and c quarks. The 
top quark is an interesting case, because in GenEvA, the top would both decay via t — > Wb 
but also radiate through t —> tg. Both of these processes need to coexist in GenEvA in order 
to cover all of multiplicity, flavor, and phase space, and this dual nature of the top may be 
interesting for understanding Monte Carlo effects in top-quark mass measurements. 



C. Choice of Evolution Variable 

For the GenEvA algorithm to be relevant in the LHC era, it will have to have a smooth in- 
terface with a hadronization/underlying event model. Regardless of the phase space strategy, 
every parton-based Monte Carlo program eventually has to interface to some external pro- 
gram that is responsible for hadronizing the initial and final state partons into the hadrons 
observed in the experiment and also to deal with secondary interactions, beam remnants, 
pile up, and other nonperturbative effects present in any hadronic collision. 



In Fig. [LJ we argued that because GenEvA is based on a parton shower and because there 
already exists a smooth interface between the showering regime and the hadronization regime 
in existing programs, then GenEvA should inherit the same smooth interface. However, as 
mentioned already, the internal GenEvA shower is optimized for speed and is not really 
meant to provide the best possible modeling of small angle QCD emissions. For example, 
the GenEvA algorithm currently evaluates a s at a fixed scale and does not implement color 
coherence. Also, it uses virtuality as the evolution variable and not p± as modern showers 
do. While it is in principle possible to construct an analytic parton shower which includes 
all these effects, this might not be the most efficient strategy. 



As discussed in Sec. |3 EL instead of using the internal shower in the showering regime, 



GenEvA can rely on an external parton shower as long as the external parton shower starts 
evolving from the event-defined matching scale t ma tch- From a physics point of view, it is 
important that the double-logarithmic dependence on t matc h cancel between the partonic and 
showering regimes, but this is difficult if the evolution variables differ in the two regimes. 
In this paper, we have been using virtuality as the matching variable, so the current version 
of the GenEvA algorithm can therefore only be smoothly interfaced with a virtuality-ordered 
parton shower. How then, can GenEvA interface with more modern p^-ordered showers? 

This turns out to be an algorithmic as opposed to a conceptual problem, though. All 
GenEvA requires to map out phase space is a parton shower with no dead zones for which 
an analytic formula exists for calculating the probability for generating a given point in 
phase space. As a phase space generator, a p_i_-ordered shower is no better or worse than a 
virtuality-ordered shower, because at the level of phase space, p^-ordering is equivalent to 
virtuality-ordering with an additional veto to eliminate certain unwanted shower histories. 

The main problem with a p_i_-ordered shower is that there is no "local" way to deter- 
mine whether a parton shower history is pj_-ordered or not. With a virtuality-ordered 
shower, every way of combining four-momenta that is consistent with flavor gives a valid 



shower history, allowing GenEvA to use the numerically efficient ALPHA algorithm |77], |78J to 
"count" the number of parton shower histories that give identical final sate four-momenta. 
It is impossible to use the ALPHA algorithm to count p^-ordered shower histories (or even 
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virtuality-ordered shower histories with an angular-ordering veto), because whether or not 
a history is p_i_-ordered (or angular-ordered) requires information about previous branches. 

It would be straightforward to rewrite GenEvA to use a p_i_-ordered shower to map out 
phase space as long as one could find a way to efficiently count p^-ordered shower histories. 
One idea is to generate the full p^-ordered expression for ^i^K^i); an d then use some 
variant of the ' Mega algorithm to generate numerically efficient code to calculate this 
sum. In other words, we can treat Yli^i^i) as a "matrix element" by itself, and construct 
an expression for it offline. 

One additional subtlety comes from the new interleaved evolution introduced in Pythia 
||. In order for reweighting to make sense, GenEvA assumes that any external showering 
program will not affect the existing GenEvA shower history. In interleaved evolution, the pri- 
mary interaction is not independent of secondary ones and the underlying event competes 
with the hard scattering for the momentum fraction of the incoming partons. Again, this is 
an algorithmic and not a conceptual problem. One could always write a more complicated 
version of GenEvA around an interleaved shower, and as long as the probability to obtain 
a given (interleaved) shower history is known, the GenEvA algorithm will be more or less 
unchanged. Practically speaking, though, the matching scale can probably be taken suf- 
ficiently large such that Pythia interleaved evolution starting from the GenEvA cutoff will 
give a reasonable description of the hadronic regime. 



D. Outlook 

Once the technical issues involving hadronic collisions, heavy resonances, and evolution 
variables are solved, there is a final technical question about how the GenEvA algorithm 
will implement more advanced matrix elements, such as those that involve complicated loop 
calculations |82| , or those based on effective theories JF0|, [71], |72fl . As currently written, the 



GenEvA algorithm assumes that there is a unique answer for the value of <r($), but depending 
on the analytic form of cr($), it may be more convenient to use Monte Carlo methods not 
only to sample phase space but also to determine er($). While the specific implementation 
will affect the performance of the GenEvA algorithm, the idea of additional "hidden" Monte 
Carlo variables being used to determine cx(<l>) is still compatible with the GenEvA algorithm's 
use of reweighting. GenEvA can in principle implement any desired differential cross section 
as long as it is well-defined on paper. 

Combined with the GenEvA framework PJ, the GenEvA algorithm has the potential to 
significantly improve our understanding of Standard Model backgrounds. The ability to 
conceive of an event simultaneously as a shower history and as a fixed-body phase space 
point allows one to answer questions that are difficult to address in either language indi- 
vidually. The efficiency of the GenEvA algorithm makes it a worthy competitor to existing 
methods, and the automatic inclusion of leading-logarithmic effects raises the bar for what 
the experimental collaborations can expect from the Monte Carlo community in the LHC 
era. 
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